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Abstract 


A  strict  kinetic,  two-dimensional  model  of  the  electron  kinetics  within  a  glow 
discharge  positive  column  is  developed.  The  problem  is  solved  in  cylindrical  geometry 
using  the  standard  two-term  Legendre  expansion  of  the  electron  velocity  distribution 
function.  The  model  establishes  a  steady  state  solution,  such  that  the  net  ionization  rate  is 
exactly  balanced  by  the  wall  loss.  In  addition  to  a  thorough  analytic  development,  we 
present  the  numerical  techniques  used  to  solve  the  resulting  elliptic  partial  differential 
equation,  including  an  efficient  method  to  treat  sparse  banded  matrices.  The  model  is 
validated  against  published  results,  local  and  nonlocal  kinetic  approximations,  and  a 
previous  Monte  Carlo  treatment.  Having  created  a  working  model,  we  conduct  an 
investigation  into  current  flow  within  the  solution  area  of  a  neon  column,  made  possible 
by  this  2-D  treatment.  Furthermore,  we  investigate  the  range  of  applicability  of  the 
earlier  local  and  nonlocal  kinetic  approximations  and  finally  present  a  short  discussion  on 
the  effect  different  forms  of  wall  loss  have  on  the  overall  distribution  function. 
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Numerical  Solutions  to  the 


Two  Dimensional  Boltzmann  Equation 

I.  Introduction 

Based  on  its  ubiquitous  application  in  neon  and  florescent  lamps,  and  more  recent 
use  in  electric  discharge  lasers,  one  might  expect  the  underlying  glow  discharge 
phenomenon  to  be  thoroughly  understood  and  computationally  well  modeled.  We  must, 
of  course,  have  an  accurate  model  of  these  devices  in  order  to  optimize  their  design 
efficiency.  However,  until  recently,  none  of  the  available  models  accurately  depicted  the 
often  dramatic  effect  of  radial  inhomogeneities  within  the  glow  discharge  device.  This 
research  addresses  and  resolves  that  deficiency  by  implementing  the  latest  modeling 
technique  -  a  two-dimensional  solution  of  the  collisional  Boltzmann  equation  (CBE), 
which  we  call  the  strict  kinetic  solution1. 

Prior  to  the  strict  solution,  most  researchers  relied  on  two  approximate  kinetic 
methods.  These  solution  methods,  known  as  the  local  and  nonlocal  approximations, 
solve  the  CBE  in  the  limit  of  high  and  low  pressure,  respectively.  Until  the  1970s  the 
local  (high-pressure)  approximation  was  used  almost  exclusively,  even  though  the 
assumptions  were  inherently  incorrect  in  most  realistic  situations.  In  the  past  twenty 
years  the  nonlocal  approximation  gained  favor  as  a  more  realistic  model  for  low-pressure 
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plasmas.  Yet  most  of  the  experimental  and  commercial  devices  in  use  today  operate  in  a 
parametric  regime  between  the  two  where  neither  the  local  or  nonlocal  approximation  is 
valid. 

Recent  advances  in  computing  power  have  facilitated  the  use  of  a  strict  kinetic 
solution  that  solves  the  complete  CBE  in  two  dimensions.  By  avoiding  the  assumptions 
inherent  to  the  approximate  methods,  the  strict  solution  is  able  to  accurately  treat  spatial 
inhomogeneities  and  remain  valid  across  a  broad  spectrum  of  pressures.  The  strict 
solution  is  in  fact  the  most  rigorous  kinetic  treatment  available  today. 

While  the  strict  solution  is  an  invaluable  tool  to  probe  the  kinetics  of  the  glow 
discharge,  the  Air  Force  is  interested  in  it  primarily  as  a  first  step  toward  much  more  lofty 
goals.  Soviet  research,  extending  over  the  past  twenty  years  and  summarized  by  Hilbun 
[1],  has  revealed  anomalous  shock  structure  and  wave  propagation  in  the  presence  of 
weakly  ionized  gases.  These  results  recently  spurred  experimental  investigations  in  the 
U.S.  that  have  confirmed  the  Soviet  findings  [2,  3];  they  include  a  reduction  in  shock 
strength,  increased  shock  standoff  distance,  and  a  general  broadening  of  the  shock 
structure.  Understanding  and  application  of  these  phenomena  could  lead  to  a  new  class 
of  hypersonic  vehicles  that  utilize  plasma  effects  to  generate  shocks  characteristic  of  one 
half  the  true  Mach  number  with  a  correspondingly  significant  decrease  in  vehicle  drag! 
However,  prior  to  designing  a  futuristic  hypersonic  fighter,  we  must  thoroughly  explore 
the  plasma-aircraft  interaction.  This  necessitates  a  refined  plasma  modeling  capability 
that  is  able  to  solve  the  CBE  in  at  least  two  dimensions.  Furthermore,  the  model  must  be 

1  The  term  “kinetic”  refers  to  any  technique  that  solves  the  CBE  directly.  “Strict”  refers  to  solving  the 
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benchmarked  in  a  parametric  regime  where  extensive  experimental  data  are  available;  the 
glow  discharge  satisfies  this  requirement.  With  this  goal  in  mind,  the  strict  solution 
detailed  in  this  thesis  represents  a  first  generation  of  the  required  modeling  capability. 

Optimistic  extrapolations  aside,  the  main  focus  of  this  research  is  to  develop, 
implement,  and  test  a  numerical  solution  to  the  spatially  inhomogeneous  CBE  -  the  strict 
kinetic  solution.  As  such,  the  body  of  this  thesis  focuses  on  the  theoretical  and  numerical 
framework  behind  the  strict  solution,  as  well  as  the  specific  details  of  the  model.  The 
remaining  chapters  follow  the  model  through  its  background,  development, 
implementation,  and  validation. 

Chapter  II  covers  background  material  leading  to  the  strict  solution.  It 
specifically  discusses  the  physics  of  the  glow  discharge,  as  well  as  giving  a  brief 
summary  of  the  techniques  and  approximations  previously  used  to  model  it.  In  Chapter 
III,  we  develop  the  strict  solution  method;  the  first  half  covers  a  theoretical  treatment  of 
the  relevant  equations,  while  the  second  discusses  the  numerical  techniques  used  to  solve 
them.  The  model  takes  shape  in  Chapter  IV;  here  we  present  “status  checks”  at  four 
phases  of  development,  using  analytic  equations  and  physical  reasoning  to  establish  the 
solution’s  validity  and  accuracy.  The  next  chapter  is  devoted  to  a  validation  of  the 
complete  model,  comparing  our  solution  to  others.  Chapter  VI  includes  a  short 
investigation  into  some  interesting  topics  made  possible  by  the  strict  solution.  Finally, 
the  last  chapter  discusses  my  ideas  for  future  work,  considering  both  expansion  of  the 
current  model  and  modifications  to  treat  related  subjects. 

CBE  completely,  avoiding  the  local  or  nonlocal  approximations. 
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II.  Background 


Elementary  physics  of  the  direct  current  (dc)  glow  discharge2  is  the  subject  of 
numerous  textbooks  (see  Nasser  [4]  and  Howatson  [5]),  and  we  can  only  briefly  review 
the  basic  concepts  here.  The  first  section  in  this  chapter  is  an  introduction  to  the  glow 
discharge  and  the  positive  column  in  particular.  The  following  two  sections  look  more 
closely  at  previous  kinetic  approximations  used  to  model  the  positive  column  and  discuss 
their  assumptions  and  range  of  applicability. 

2.1.  Introduction  to  the  Positive  Column 

2.1.1.  Glow  Discharge  Basics 

A  simple  glow  discharge  device  consists  of  a  cylindrical  glass  tube,  filled  with  a 
gas  (see  Figure  1).  Most  tubes  are  one  to  two  centimeters  in  diameter,  and  between  ten 
and  one  hundred  centimeters  long.  Within  the  tube,  the  gas  pressure  normally  ranges 
from  a  few  hundredths  to  tens  of  a  Torr.  The  discharge  is  initiated  and  sustained  by  an 
external  dc  power  supply.  The  discharge  voltage  is  typically  several  hundred  Volts,  with 
corresponding  currents  of  tenths  to  hundreds  of  milli-amperes  [6].  The  power  is  coupled 
to  the  plasma  through  two  electrodes,  one  at  each  end  of  the  tube. 

In  very  simple  terms,  the  glow  discharge  is  created  as  the  externally  provided 
power  is  converted  into  light  and  heat.  The  conversion  occurs  when  electrons, 
accelerated  by  the  applied  axial  electric  field,  collide  with  the  neutral  gas.  Some  of  these 

2  This  study  considers  only  the  direct  current  discharge.  Suitable  modifications  to  the  model  would  allow 
us  to  account  for  an  alternating  current  as  well. 
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energetic  electrons  have  enough  energy  to  excite  the  neutrals  into  higher  electronic  states. 
These  neutrals  in  turn  emit  radiation  as  they  relax  to  the  ground  state.  This  radiation  is 
the  glow. 


Cathode  Negative  Positive  Anode  Anode 


Figure  1  Schematic  of  a  typical  neon  discharge  in  a  50  cm  tube  at  1  Torr.  Luminous  regions  are 
shown  shaded,  along  with  plots  of  luminosity,  electric  field,  and  charge  density  [4,  pg  399]. 

Figure  1  illustrates  a  schematic  of  a  typical  glow  discharge  device,  as  well  as 
relative  values  of  important  parameters  along  the  length  of  the  tube.  Studying  the  figure, 
we  note  that  the  glow  discharge  can  be  divided  into  a  number  of  physically  meaningful 
regions  [4,  pp  397-425].  At  the  cathode  end,  the  impact  of  heavy  ions  onto  the  cathode 
surface  releases  electrons  into  the  system.  At  the  time  of  their  release,  the  electrons  do 
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not  have  enough  energy  to  excite  electronic  levels  in  the  neutrals;  hence  there  is  no 
radiative  decay  and  no  glow  emission  in  this  region,  known  as  the  Aston  Dark  Space. 

The  strong  axial  electric  field  accelerates  the  freed  electrons  toward  the  anode, 
and  their  kinetic  energy  soon  exceeds  the  excitation  threshold  of  the  neutral  species.  This 
marks  the  beginning  of  the  Cathode  Glow.  Note  that  in  this  region  the  luminosity  is  still 
quite  low,  since  the  number  of  electrons  is  relatively  low.  The  electrons  continue  to 
accelerate,  achieving  energies  in  excess  of  the  neutral’s  ionization  threshold  at  the 
Cathode  Dark  Space  boundary.  When  this  occurs,  ionization  begins  to  compete  with 
excitation  for  electron  energy,  reducing  the  overall  luminosity.  The  luminosity  does  not 
go  to  zero  however;  in  fact  all  the  dark  regions  have  some  luminosity,  just  at  a  lower 
level  than  their  neighbors. 

The  electrons  released  by  ionization  in  the  Cathode  Dark  Space  are  accelerated  to 
excitation  energies,  and  because  of  the  increase  in  electron  density,  the  Negative  Glow  is 
the  brightest  band  in  the  entire  tube.  The  large  number  of  electrons  in  this  region  causes 
a  weakening  in  the  axial  field  strength,  and  it  can  actually  reverse  for  a  short  distance  [4, 
pg  412].  Because  of  the  weak  field,  the  acceleration  in  this  region  is  not  strong  enough  to 
keep  the  electrons  above  the  excitation  threshold,  and  the  Faraday  Dark  Space  results.  In 
the  Faraday  Space,  the  number  of  electrons  tapers  off,  and  the  axial  electric  field 
increases,  until  it  is  once  again  able  to  accelerate  electrons  to  excitation  energies. 

At  this  point,  the  Positive  Column  begins,  maintaining  a  nearly  uniform  glow. 
Near  the  end  of  this  column,  the  axial  field  begins  to  increase,  and  once  again  the 
electrons  gain  enough  energy  for  ionization  to  strongly  compete  with  excitation.  The 
final  two  regions  are  analogous  to  the  Cathode  Dark  Space  and  Negative  Glow.  In  the 


6 


Anode  Dark  Space  the  luminosity  decreases  due  to  ionization;  in  the  Anode  Glow  newly 
formed  electrons,  accelerated  by  the  increasing  electric  field,  again  reach  the  excitation 
threshold. 

Returning  to  the  positive  column  region  in  Figure  1,  we  note  a  number  of 
interesting  features;  the  axial  electric  field,  luminosity,  electron  number  density,  and 
current  are  all  essentially  constant.  In  addition,  while  the  size  of  the  other  regions 
depends  on  the  electron  mean  free  path,  the  length  of  the  positive  column  scales  with  the 
tube.  As  the  tube  grows  longer,  so  does  the  positive  column;  if  the  tube  is  too  short,  the 
column  may  disappear  altogether.  The  reason  for  this  scaling  comes  from  its  purpose  -  it 
provides  electrical  continuity  between  the  cathode  and  anode  [6].  Because  the  positive 
column  grows  with  tube  length,  it  is  often  the  largest  of  the  glow  discharge  segments. 
These  properties  will  be  exploited  later  to  simplify  our  model  of  the  glow  discharge. 

Although  Figure  1  and  the  resulting  discussion  considered  only  axial 
characteristics  of  the  discharge  tube,  radial  variations  exist  as  well.  In  Figure  2,  we  see 
plots  of  experimentally  determined  radial  space  charge  potentials  within  the  cylindrically 
symmetric  positive  column.  These  potentials  form  as  relatively  light  electrons  speed 
toward  the  tube  wall,  leaving  the  heavier  ions  behind.  This  charge  separation  creates  a 
potential  and  hence  a  radial  component  to  the  electric  field.  Because  of  the  radial 
potential,  any  model  which  hopes  to  accurately  describe  the  glow  discharge  needs  to 
consider  not  only  axial  variations,  but  radial  ones  as  well. 
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Figure  2  Experimentally  determined  radial  potential  (— V[r]  in  Volts)  as  a  function 
of  normalized  radius  for  neon  and  krypton  at  .62  and  4.8  Torr  respectively  [7], 

Due  to  the  shear  numbers  of  particles  involved,  electrons  in  the  system  are  often 
treated  statistically  and  characterized  using  a  distribution  function.  The  electron  velocity 
distribution  function  (EVDF),  labeled/,  represents  the  number  of  particles  found  at  a  set 
of  coordinates,  e.g .f[x,  y,  z,  vx,  vy,  vz,  t].  To  find  the  total  number  of  particles  at  a  given 

spatial  location  the  EVDF  is  integrated  over  velocity  space,  i.e.  f°  f[r,v]dv  =  n  [r] . 

J—oo  e 

By  considering  an  electron  distribution,  we  ignore  the  motions  of  individual 
particles,  treating  the  plasma  statistically.  This  approximation  is  justified  for  all  but  the 
most  tenuous  of  gases,  and  it  is  certainly  acceptable  in  our  pressure  range  of  interest,  0.1- 
10  Torr.  Once  the  distribution  function  is  known,  we  can  derive  all  of  the  relevant 
macroscopic  quantities  using  moments  of  the  Boltzmann  equation  (see  Appendix  A). 
These  quantities,  such  as  the  excitation  collision  rate  (equation  (44)),  correspond  directly 
to  important  measurable  parameters  like  the  luminosity. 
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Often  the  distribution  function  is  transformed  from  velocity  to  energy  coordinates, 
and  then  referred  to  as  the  electron  energy  distribution  function  (EEDF).  Figure  3  shows 
examples  of  four  different  EEDFs.  The  first  two,  (a)  and  (b),  result  from  approximate 
solutions  of  the  CBE,  under  the  assumption  of  constant  collision  frequency  and  constant 
collision  cross-section  respectively.  The  third  and  fourth  represent  more  physically 
realistic  distributions,  for  the  case  of  a  noble  gas,  (c),  and  a  simple  molecular  gas,  (d). 


Figure  3  Four  examples  of  the  electron  energy  distribution  function  plotted  on  a  log/linear 
scale:  (a)  Maxwellian,  (b)  Dmyvesteyn,  (c)  typical  Noble  Gas,  (d)  typical  Molecular  gas. 

2.1.2.  Modeling  the  Glow  Discharge 

When  trying  to  realistically  model  the  electron  kinetics  of  a  glow  discharge,  we 
are  forced  to  make  a  number  of  significant  simplifications.  The  first  simplification  is  to 
recognize  that  while  most  of  the  glow’s  segments  are  fixed  in  length  for  a  given  set  of 
conditions,  the  positive  column  grows  and  shrinks  with  the  tube  length  (section  2.1.1).  If 
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the  discharge  tube  is  more  than  a  few  centimeters  long,  the  positive  column  accounts  for 
the  majority  of  the  discharge,  and  is  thus  considered  the  working  section  for  most 
applications.  To  simplify  our  study,  this  treatment  considers  only  the  positive  column3. 

Having  chosen  to  model  only  the  positive  column  region  of  the  glow  discharge, 
we  can  then  make  use  of  its  cylindrical  symmetry  and  axial  homogeneity  (section  2.1.1). 
Accounting  for  these  two  properties,  the  required  spatial  dimension  of  the  solution  is 
reduced  from  three  to  one  -  the  radial  dimension. 

Even  with  this  simplification,  modeling  the  dc  positive  column  is  extremely 
complicated;  past  authors  have  used  a  number  of  different  techniques  to  tackle  the 
problem.  These  techniques  fall  into  three  broad  classes:  Moment,  Particle,  and  Kinetic 
methods.  The  moment  methods  include  such  classic  treatments  as  Schottky’s  ambipolar 
diffusion  theory  and  Tonks  and  Langmuir’s  free-fall  theory.  As  the  name  suggests, 
moment  methods  use  moments  of  the  Boltzmann  equation,  treating  the  electrons  as  a 
fluid.  These  methods  start  with  an  assumed  form  for  the  EVDF,  and  use  the  resulting 
transport  coefficients  and  collision  frequencies  to  transform  the  problem  into  a  set  of 
fluid  balance  equations.  Yet  without  experimental  verification  of  the  assumed 
distributions,  the  results  are  only  qualitative.  In  spite  of  this,  historically,  most  theoretical 
treatments  of  the  positive  column  have  relied  on  the  moment  method  [8]. 

Particle  methods,  which  include  Monte  Carlo  and  Particle-in-Cell,  consider  the 
forces  and  collisions  acting  on  individual  particles  to  predict  their  motions.  Since  it  is 
currently  impossible  to  account  for  the  movement  of  over  1016  particles/cm3  (~  1  Torr), 

3  Cathode  to  anode  models  exist  which  attempt  to  account  for  the  physics  along  the  entire  length  of  the 
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these  techniques  attempt  to  create  a  representative  sample,  and  apply  probability  theory 
to  extrapolate  the  total  answer.  Even  so,  particle  methods  require  significant 
computational  effort  to  obtain  a  solution,  and  even  then  are  subject  to  sampling  errors. 
Due  to  the  computation  time  required,  recent  authors  have  suggested  that  techniques  such 
as  Monte  Carlo  are  best  suited  to  serve  as  benchmarks  for  other,  faster  methods  [9]. 

The  third  class,  the  kinetic  method,  is  the  subject  of  this  thesis.  Kinetic  methods 
solve  directly  for  the  EVDF  using  the  collisional  Boltzmann  equation, 

+  f  +  a-V /  =  (— )  ^ 

where  /  =  f[r,v,t]  is  the  EVDF. 

The  Boltzmann  equation  is  in  essence  a  continuity  equation  in  both  configuration 
and  velocity  space.  The  first  term  in  equation  (1)  represents  temporal  changes  to  the 
EVDF,  which  for  steady  state  solutions  is  zero.  The  second  and  third  terms  denote 

convection  in  configuration  and  velocity  space  respectively,  where  Vr  and  Vv  are  the 

corresponding  gradient  operators.  The  acceleration,  a ,  contains  all  of  the  forces  acting 
on  the  electron,  which  in  our  study  includes  both  a  radial  and  axial  electric  field.  The 

final  term,  (J-)collision  ,  accounts  for  all  of  the  collisions.  It  acts  as  both  a  source  and  a 
at 

loss,  depending  on  the  type  of  collision  and  region  of  velocity  space  considered.  A 
discussion  of  the  various  collisions  included  in  this  model  is  presented  in  section  3.1.2. 

Because  of  the  difficulty  associated  with  solving  equation  (1)  directly,  early 
pioneers  of  the  kinetic  method  looked  to  further  simplify  the  problem.  They  did  this  by 

discharge  tube,  but  this  treatment  considers  only  the  positive  column. 
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considering  limiting  forms  for  the  CBE  in  the  high  and  low-pressure  range.  These  two 
limits  correspond  respectively  to  the  local  and  nonlocal  kinetic  approximations.  While 
the  treatment  in  this  thesis  uses  a  strict  kinetic  solution,  which  is  valid  in  both  pressure 
regimes,  we  will  gain  significant  insight  by  first  reviewing  the  approximate  solutions. 

2.2.  The  Local  Approximation 


One  of  the  first  solutions  using  the  kinetic  method  was  published  in  1946  by 
Holstein  [10].  In  this  work,  Holstein  investigated  the  properties  of  a  discharge  in  which 
he  assumed  the  electric  field  to  be  everywhere  homogeneous.  Under  these  conditions, 
the  spatial  gradient  of  the  distribution  function  must  be  zero,  and  the  expansion  of 
equation  (1)  results  in  an  ordinary  differential  equation  (ODE).  The  ODE  is  given  below, 


d  MeJ^f  d 
d£ K  3 NQ[e]  d£ 


F[£])  + 


^-^-(£2NQ[£]F[£])  +  S'  =  0 
Mn  d£ 


(2) 


where  N  is  the  neutral  number  density,  Q[e\  the  total  collision  cross  section,  and  S' 
accounts  for  inelastic  collisions.  F[e]  is  the  EEDF,  normalized  such  that 


^F[£]£V2d£  =  1. 


If  we  consider  only  certain  types  of  collisions,  equation  (2)  can  be  solved 
analytically.  Neglecting  inelastic  collisions  and  assuming  a  constant  elastic  collision 
frequency,  the  analytic  result  is  the  Maxwellian  distribution  shown  in  Figure  3  a.  If,  on 
the  other  hand,  elastic  collisions  are  represented  with  a  constant  collisional  cross  section, 
the  solution  is  a  Druyvesteyn  distribution  given  by, 
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where  Qel  refers  now  to  the  elastic  cross  section  only.  The  Druyvesteyn  distribution  is 
shown  in  Figure  3b.  Although  these  analytic  solutions  exist,  in  practice  equation  (2)  is 
solved  numerically  using  a  complete  set  of  experimentally  consistent  elastic  and  inelastic 
cross-sections  [11]. 

While  the  assumption  of  a  homogeneous  electric  field  greatly  simplifies  the  CBE, 
it  is  rarely  true  in  a  physical  discharge  device.  As  discussed  in  section  2.1.1,  we  can  treat 
the  positive  column  as  having  a  homogeneous  axial  electric  field,  but  Figure  2 
demonstrates  that  radial  inhomogeneities  do  exist.  However,  under  certain  conditions,  it 
is  possible  to  neglect  the  effect  of  these  radial  variations,  and  justify  the  use  of  equation 
(2).  If  the  pressure  in  the  discharge  is  high  enough,  the  electrons  will  undergo  sufficient 
collisions  to  lose  all  of  their  energy  (relax  in  energy  space)  prior  to  traveling  an 
appreciable  distance  radially.  Assuming  the  radial  distance  traveled  is  small  enough,  the 
electrons  will  not  “see”  any  variation  in  the  electric  field.  Under  these  conditions,  the 
energy  distribution  of  electrons  is  determined  solely  by  the  local  electric  field  and 
equation  (2)  is  justified.  This  assumption  is  called  the  Local  Approximation. 

In  order  for  the  local  approximation  to  be  accurate,  the  electron  energy  relaxation 
length  must  be  much  shorter  than  the  characteristic  length  for  field  variations.  One  way 
to  quantify  this  requirement  is  through  the  PR  value,  where  P  is  the  pressure  in  Torr,  and 
R  is  the  tube  radius  in  centimeters4.  Using  a  random  walk  analysis  and  the  fact  that  an 


4  A  more  appropriate  measure  would  be  NR  where  N  is  the  number  of  molecules  and  R  the  radius,  but  I 
will  follow  the  convention  of  Ingold  [8]  and  use  PR. 
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electron  loses  only  a  small  fraction, - ,  of  its  energy  per  elastic  collision,  the  electron 

Mn 

energy  relaxation  length  L  is  given  by 

L2\=MjLN-2Q-2  (4) 

/  r\  x-' momentum 

Zm  transfer 

where  Q  is  the  average  momentum  transfer  cross-section  in  angstroms  squared  and  N  the 
neutral  number  density.  Taking  the  characteristic  distance  for  field  variations  to  be  the 
tube  radius  R,  we  can  approximate  the  “PR”  delimiter  as 

NR  »  60 /T  yjMjm  Q~fomentum 

transfer 

PR  »  0.2 N/m  Qmomentum  ^ 

transfer 

where  T  is  300K.  Applying  equation  (5)  to  neon  gas,  the  relaxation  length  requirement  is 
met  if  PR  »  17  Torr-cm.  Given  that  lamps  and  electric  discharge  lasers  operate  at  PR 
values  on  the  order  of  one  Torr-cm,  this  requirement  is  rarely  satisfied  in  practice. 

Further  difficulties  exist  in  applying  the  local  approximation  to  realistic  devices. 
In  deriving  equation  (5)  we  chose  the  characteristic  distance  for  field  variations  to  be  the 
tube  radius.  But  as  shown  in  Figure  2,  the  radial  potential,  and  hence  its  derivative,  the 
radial  electric  field,  vary  rapidly  near  the  tube  wall.  These  rapid  changes  invalidate  the 
earlier  approximation  of  the  electrons  not  “seeing”  field  variations.  We  will  find  that 
even  when  the  requirement  of  equation  (5)  is  satisfied,  the  local  approximation  is  never 
accurate  near  the  tube  wall. 

2.3.  The  Nonlocal  Approximation 

When  the  PR  value  is  much  less  than  equation  (5),  the  electrons  are  able  to 
repeatedly  traverse  the  discharge  tube  before  relaxing  in  energy  space.  For  this  low- 
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pressure  regime,  the  electrons  sample  the  entire  tube  and  carry  with  them  a  “memory”  of 
the  electric  field  variations  throughout.  Because  of  this,  their  energy  distribution  is  no 
longer  a  function  of  the  local  electric  field,  rather  it  is  a  function  of  the  average  electric 
field  and  its  variation.  In  this  case  a  new,  approximate  solution  is  available  -  the 
Nonlocal  Approximation.  Bernstein  and  Holstein  were  the  first  to  explore  this  low- 
pressure  limit,  publishing  a  paper  in  1954  [12]. 

Taking  into  account  the  radial  potential  shown  in  Figure  2,  the  expansion  of 
equation  (1)  results  in  a  partial  differential  equation  (PDE)  in  two  dimensions,  energy  and 
radius.  Much  like  the  local  case,  the  nonlocal  approximation  simplifies  the  problem  to  an 
ODE  in  energy  only.  In  this  instance,  however,  the  transformation  is  achieved  by 
averaging  over  the  spatial  derivatives,  rather  than  neglecting  them  [6].  In  order  to  justify 
the  radial  averages,  the  electron  energy  relaxation  length  must  greatly  exceed  the  tube 
radius.  In  other  words,  the  pressure  or  PR  of  the  system  must  be  low,  which  is  the 
opposite  of  the  condition  in  equation  (5). 

After  averaging,  the  CBE  becomes 


j~{  ^^~(P[e]F[e])  V  (<?[£]  ^[e])  +  S'  =  0 

de  3 N  de  MN  de 


p.fe-fr))  rdr 

Jo  Q[e-(p[r]] 

q[e]  =  £  (e  -  (p[r)f  Q[£  -  (p[r]y  dr 


Q[£-(p[r]] 


where  Q  is  again  the  total  collision  cross  section,  but  S"  is  now  an  averaged  inelastic 
collision  term.  The  two  additional  terms,  p[e ]  and  q[e],  calculate  the  necessary  averages 
over  the  dimensions  of  the  tube,  and  involve  the  electron  motive,  (p[r]  =  -V[r] .  With  the 


addition  of  a  radial  potential,  the  kinetic  energy,  u,  depends  not  only  on  e,  but  r  as  well. 
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Equation  (6)  looks  very  similar  to  (2)  for  the  local  approximation,  and  as  in  that 
case  analytic  solutions  exist  under  certain  conditions.  Assuming  no  inelastic  collisions 
and  a  constant  elastic  cross  section  we  arrive  at  a  modified  Druyvesteyn  distribution, 


Nonlocal  l 
Dru 


T[3/4]  M 


-2m(NQel)2£2 

MN(e0Ezf 


where  we  have  assumed  a  quadratic  potential  (the  solution  is  slightly  different  depending 
on  the  form  of  the  potential).  Bernstein  and  Holstein  [12]  first  derived  equation  (7)  in 
planar  geometry,  but  we  have  converted  it  to  cylindrical  coordinates  for  this  treatment. 

While  the  local  and  nonlocal  approximations  adequately  describe  discharge 
kinetics  in  their  respective  pressure  regimes,  the  regimes  are  not  well  defined.  In 
addition,  there  exists  a  significant  gap  between  the  two  approximations  where  neither  is 
very  accurate.  Most  physical  glow  discharge  devices  live  in  this  gap  and  are  thus,  at 
present,  poorly  described.  Only  through  a  strict  solution  of  the  complete,  expanded  CBE 
can  we  accurately  model  the  middle  pressure  regime.  The  remainder  of  this  thesis 
investigates  that  strict  kinetic  solution. 
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III.  Methodology  of  the  Strict  Kinetic  Solution 


This  chapter  develops  the  theoretical  equations  used  to  model  the  positive  column 
region  of  a  dc  glow  discharge,  then  presents  the  numerical  techniques  used  to  solve  the 
resulting  partial  differential  equation.  The  development  in  section  3.1.1  is  necessarily 
condensed,  and  focuses  on  understanding  the  assumptions  made  and  resulting  terms  in 
the  equations.  A  more  thorough  derivation  is  provided  in  Appendix  B. 

3.1.  Theoretical  Development  of  Solution  Equations 

3.1.1.  Expansion  of  the  Collisional  Boltzmann  Equation 

Using  the  strict  kinetic  method,  we  solve  for  the  distribution  of  electrons  within 
the  positive  column  as  a  function  of  position  and  velocity  (or  energy).  As  stated  in 
section  2.1.2,  the  starting  point  for  the  derivation  is  the  collisional  Boltzmann  equation, 

|-+V'Vr/  +  a.V,/  =  (|-)-a„1,  (8) 

where  /  =  f[r ,  v,t]  is  the  electron  velocity  distribution  function  normalized  to  the 

electron  density,  i.e.  j~  f[r,v]dv  =  ne[r] . 

Reiterating  the  discussion  of  section  2.1.2,  the  Boltzmann  equation  is  an  electron 
continuity  equation  in  both  configuration  and  velocity  space.  In  our  steady-state 
treatment,  the  first  term,  representing  temporal  changes  to  the  EVDF,  is  zero.  The 
remaining  terms  balance  the  divergence  of  the  electron  flux  with  gains  and  losses  due  to 
collisions.  A  complete  discussion  of  the  various  collisions  included  in  the  model  is 
presented  in  section  3.1.2. 
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By  writing  a  Boltzmann  equation  only  for  the  electrons  and  not  the  other 
“secondary”  constituents  (e.g.  ground  state  neutrals,  ions,  etc.),  we  are  essentially 
assuming  a  form  for  their  distributions.  In  this  treatment  we  treat  the  secondary  species 
as  being  uniformly  distributed  spatially,  and  as  having  negligible  temperature.  This  is  a 
reasonable  assumption  given  the  large  mass  ratio  between  the  secondary  constituents  and 
electrons,  and  considering  the  relatively  low  neutral  gas  temperatures  observed  within  the 
plasma.  Making  it  allows  us  to  de-couple  an  otherwise  complicated  system.  We  are  left 
with  the  single  Boltzmann  equation,  albeit  in  six  dimensions,  describing  an  EVDF  that 
contains  all  of  the  phenomenon  of  interest  in  the  system  [13]. 

Solving  the  Boltzmann  equation  in  six  dimensions  is  a  daunting  task,  so  further 
assumptions  are  made  to  simplify  the  problem.  The  next  widely  used  assumption  is  to 
assert  that  for  the  conditions  of  interest  the  EVDF  is  nearly  isotropic.  This  allows  us  to 
expand  the  EVDF  in  spherical  harmonics  in  velocity  space,  keeping  only  the  first  two 
terms. 


/[/in  =  XF,(cos[0])/,(r,v) 


=  f0[r,v]  + 


v-/,[r,v] 


(9) 


v  •  f. 

Here  f0  is  the  isotropic  part  of  the  EVDF,  while - L  refers  to  the  anisotropic  part. 

v 

Notice  that  /,  is  a  vector,  which  can  be  resolved  into  radial  and  axial  components,/,,  and 
fz  respectively.  Further  simplifications  are  made  by  moving  from  v  =>  v  ,  incorporating 
the  vector  information  directly  into  the  equations.  By  recognizing  the  cylindrical 
symmetry  and  axial  homogeneity  of  the  positive  column  (section  2.1.1)  we  finally  reduce 
the  problem  from  six  dimensions  down  to  two,  one  radial  and  one  speed. 
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The  previous  two-term  expansion  is  valid  if:  the  energy  gained  from  the  electric 
field  between  collisions  is  small,  and  the  total  number  of  elastic  collisions  is  significantly 
greater  than  the  inelastic  ones.  These  conditions  ensure  the  electrons  are  well  distributed 
in  velocity  space,  i.e.  isotropic,  and  for  the  conditions  of  interest  are  generally  considered 
valid  [14,  pg  2], 

As  derived  in  Appendix  B,  the  solution  continues  by  substituting  (9)  into  (8), 

expanding  the  V  operators,  and  keeping  just  the  zeroth  and  first  order  terms.  We  also 
consider  the  source/loss  terms  in  the  collision  integral,  adding  the  contributions  from 
elastic  and  inelastic  collisions,  as  well  as  electron  production  arising  from  collisions 
between  excited  state  neutrals. 

With  further  manipulations,  the  final  result  is  a  single  elliptic  PDE.  In  this 
transformed  state,  the  solution  is  a  function  of  radial  position  and  total  energy;  in  other 

words,  /0  is  now  the  electron  energy  distribution  function.  The  CBE  for  the  strict  kinetic 
method  becomes 


1  d  (  ru  d  7  ,  |  d  u(e0Ez)2  d 
r  dr  3  K[r,e]dr  0  de  3  K[r,e]  de 


fo)  +  -^(G[r,£]f0)-uH[r,e]f0  +  S0[r,e,f0]  =  0 


(10) 


where 


G[r,e]  =  Y2^u2Nk  Qekl[u]  (11) 

*  Mk 

+  +  (12) 

k  l  k  l  k 

K[r,e]  =  ^NkQekl[u]+H[r,e ]  (13) 

k 

S0[r,u,f0]  Ukl  Nk  <2r;"X,]  foir,u'u]  (14) 

k  l 
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These  term  definitions  follow  the  formalism  of  Uhrlandt  and  Winkler  [7,  pg  520]. 
Specific  parameters  include:  Ez  -  the  axial  electric  field;  Nk  -  the  number  density,  where 
k  refers  not  only  to  different  neutral  species,  but  also  to  different  excited  states  of  the 
same  species;  m  -  mass  of  the  electron;  M k  -  mass  of  the  neutral  species  k;  and 

QP™cess[u]  .  ^e  collisional  cross  section,  which  depends  not  only  on  the  process 
involved,  but  the  kinetic  energy  and  collision  partners. 

Interpreting  equations  (1 1)-(17),  G[r,e]  is  the  scaled  elastic  momentum  transfer 
term.  H[r,e]  sums  the  three  inelastic  collisions  (excitation,  de-excitation,  and  ionization), 
becoming  the  total  inelastic  term,  while  K[r,e]  includes  both  elastic  and  inelastic 
collisions,  giving  the  total  momentum  transfer  term.  S0[r,e,f0]  is  a  source  term,  adding 

electrons  at  energy  u  based  on  inelastic  collisions  taking  place  at  other  energies  u  . 

Upon  closer  inspection,  we  recognize  that  (10),  like  (8),  is  a  continuity  equation. 

The  first  term,  1— ( — — — —  7) ,  represents  the  divergence  of  the  radial  flux;  as  u 
rdr  3K[r,e]  dr  Jo 


increases,  the  flux  increases,  while  as  K[r,e]  (the  total  collision  term)  increases,  the  radial 

flux  decreases.  The  second  term,  — —  7\ ,  denotes  the  divergence  in  energy 

de  3K[r,e]  3e  0 


space  of  an  electric  field  driven  flux;  again  this  flux  tends  to  increase  as  u  increases,  yet 


diminishes  as  K[r,e ]  increases.  The  next  term,  —  (G[r,£]/0) ,  is  also  the  divergence  of  an 

de 

energy  flux,  but  this  time  due  to  elastic  collisions.  Finally,  the  last  two  terms  are  the 
losses  and  sources  due  to  inelastic  collisions;  —  uH[r,e]f0  accounts  for  the  losses  at  a 
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given  r  and  £  that  occur  due  to  the  various  inelastic  collisions.  Corresponding  gains  from 
other  energies  appear  in  S0[r,e,f0] . 

In  the  course  of  deriving  equation  (10),  the  following  relationships  for  the 


The  next  section  addresses  these  issues. 
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Figure  4  Energy  level  diagram  depicting  three  types  of  inelastic  electron 
collisions  involving  neutral  species  k  and  excited  states  /  and 

3.1.2.  Discussion  of  Collision  Terms 

Careful  consideration  of  the  relevant  collision  mechanisms  is  required  to  achieve 
a  robust  and  flexible  simulation,  yet  it  is  not  necessary  to  include  every  possible  type  of 
collision  in  order  to  create  a  quantitatively  accurate  model  of  the  system.  Many  of  the 
collisions,  such  as  electron/electron,  rarely  occur  in  the  weakly  ionized  gases  considered 
and  therefore  contribute  little  to  the  overall  electron  distribution.  By  disregarding 
negligible  terms,  the  model  remains  accurate  and  relatively  simple  to  implement. 

Our  treatment  considers  three  classes  of  inelastic  collisions:  excitation,  de¬ 
excitation,  and  ionization.  We  ignore  volumetric  recombination,  choosing  instead  to 
assume  that  all  recombination  occurs  at  the  wall  as  part  of  the  wall  loss  boundary 
condition.  This  is  a  reasonable  assumption  for  the  low  electron/ion  densities  considered 
as  well  as  the  relatively  high  electron  temperature  that  results.  In  the  case  of  elastic 
collisions,  we  include  only  electron/neutral  and  disregard  electron/ion  and 
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electron/electron  collisions.  This  again  limits  the  treatment  to  weakly  ionized  plasmas  in 
which  these  types  of  collisions  are  negligible. 


Collision  of  1st  Kind  (Excitation) 
e  +  N— >  e'  +  N* 


Ionization  Collision 
e  +  N  — »  2e  +  N+ 


u 


it  —  2  u  +  u 


ionization 

threshold 


Collision  of  2nd  Kind  (De-excitation) 
e  +  N*  — >  e'  +  N 


Collision  of  1st  and  2nd  kind  involve  shifting  one 
electron.  Ionization  creates  an  electron  and  shifts 
another. 


Figure  5  Electron  scattering  in  three  types  of  inelastic  electron  collisions. 

The  inelastic  collisions  considered  in  this  model  are  represented  in  Figure  4  and 
Figure  5.  Figure  4  illustrates  the  effect  inelastic  collisions  have  on  the  neutral  species’ 
internal  energy  levels.  The  complementary  view  in  Figure  5  shows  the  effect  of  the  same 
collisions  on  the  kinetic  energy  of  the  electrons. 

Excitation,  or  a  collision  of  the  first  kind,  inelastically  scatters  energetic  electrons 
to  lower  energies,  while  promoting  the  colliding  atom  or  molecule  to  a  higher  internal 
energy  state.  The  electron’s  energy  loss  is  equal  to  the  excitation  threshold  of  the  neutral 
species.  Because  the  plasma  in  this  treatment  is  weakly  ionized,  excitation  of  ions  is 
neglected.  The  density  of  excited  state  neutrals  is  also  relatively  small;  thus  any  further 
excitation  to  higher  electronic  levels  is  also  neglected. 

The  reverse  process  is  de-excitation,  a  collision  of  the  second  kind,  or  a  super¬ 
elastic  collision.  In  this  case,  collisions  occur  between  electrons  and  excited  state 
neutrals,  scattering  the  electrons  to  higher  energies  as  the  excitation  energy  of  the  neutral 
is  transferred  to  kinetic  energy  in  the  electron.  Once  again  the  ion  density  is  considered 
too  small  for  excited-state  ion/electron  collisions  to  contribute  to  the  solution. 
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In  the  case  of  ionizing  collisions,  we  consider  two  types:  between  energetic 
electrons  and  ground  state  neutrals  (Figure  4  and  Figure  5),  and  between  excited  state 
neutrals  (Penning  ionization).  As  in  the  case  of  excitation,  collisions  between  electrons 
and  excited  state  neutrals  (step-wise  ionization)  are  disregarded. 

The  Penning  ionization  term  is  very  important,  from  both  a  physical  and 
computational  standpoint.  Consider  the  homogenous  nature  of  equation  (10),  where 
every  term  depends  on/0;  this  means  there  are  an  infinite  number  of  proportional 

solutions.  In  contrast,  the  Penning  term  depends  only  on  the  number  density  of  excited 
state  neutrals  and  the  rate  coefficient  of  electron  production.  Introduction  of  this  term 
adds  an  inhomogeneity  to  the  equation  set,  and  resolves  the  solution.  An  inhomogeneous 
term  can  also  arise  through  the  introduction  of  certain  boundary  conditions,  but  not  all 
boundary  conditions  produce  this  effect.  It  is  therefore  important  to  include  the  Penning 
ionization  in  order  to  obtain  an  absolute  solution. 

Returning  to  equation  (10),  the  expression  for  S0  is  now  expanded  to  show  the 
contribution  of  each  inelastic  collision: 

S»[r.«./0]=ESJ?  (17) 

k  k*  V  l 

+  2X  (u+u%)  Nk  j2“[«  +  «“]  /0[r,«  +  ii“] 

k  l 

+  £4  (2 «  +  <)  Nk  Qi° [2u  +  u[°]  f0[r,2u  +  u‘k°] 

k 

+  XX  (»-<) N„  QHu-uU  /„[/%«-<] 

k  l 

Ignoring  for  a  moment  the  Penning  term,  the  other  three  are  strikingly  similar,  combining 
a  collisional  interaction  energy,  neutral  number  density,  cross  section,  and/0.  The  major 

difference  between  them  is  the  interaction  energy.  This  energy  is  the  amount  required  for 
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the  scattered  electron  to  end  up  with  kinetic  energy  u.  For  example,  under  excitation,  the 
electron  collides  with  a  neutral,  transferring  energy  equal  to  the  excitation  threshold  (see 
Figure  4  and  Figure  5).  Therefore,  for  the  electron  to  end  up  with  kinetic  energy  u  after 
the  collision,  it  must  have  energy  u  +  uexc“anon  before  the  collision.  The  rationale  is 
reversed  for  de-excitation.  In  the  special  case  of  ionization  a  new  electron  is  formed,  thus 
the  interaction  energy  must  be  sufficient  to  overcome  the  ionization  threshold  and 
provide  kinetic  energy  u  to  both  electrons.  By  giving  both  electrons  kinetic  energy  u,  we 
have  assumed  that  they  equally  share  the  excess  energy. 

Unlike  excitation  and  de-excitation,  the  ionization  term  contains  a  prefactor  of  4. 
This  4  is  actually  comprised  of  two  parts:  a  factor  of  two  to  account  for  the  newly 
released  electron  and  a  second  factor  of  two  that  relates  to  the  energy  bin  spacing. 
Consider  the  case  of  excitation,  in  that  collision  the  scattered  electron  maps  directly  from 
one  energy  bin  ( u  +  u exc,tanon )  to  another  (w).  In  ionization,  however,  the  electron  starts  in 
an  energy  bin  ( 2u  +  uwn,m“on )  that  is  twice  as  wide  as  the  destination  bin  (w).  The  second 
factor  of  two  accounts  for  this  difference  in  bin  widths. 

Returning  to  the  case  of  Penning  ionization,  the  number  of  electrons  produced  at 
energy  u  is  a  function  of  the  number  density  of  excited  state  neutrals  and  rate  coefficient 
of  electron  production,  zkk* .  The  delta  function,  8[uPenmng  ] ,  represents  the  assumption 

that  all  Penning  electrons  are  produced  at  a  single  kinetic  energy5,  uPennins.  When  the 
excited  state  neutrals  collide,  they  relax  to  their  ground  states;  the  released  energy  frees 

5  Some  authors  [7]  choose  to  spread  the  Penning  electrons  out  across  a  range  of  energies,  but  that  is  not 
done  in  this  model. 
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an  electron,  with  any  remaining  energy  going  into  electron  kinetic  energy.  Thus  uPennins 
is  equal  to  the  sum  of  excitation  thresholds  of  the  colliding  neutrals,  minus  the  ionization 
threshold  of  the  gas  species. 

3.1.3.  Boundary  Conditions  to  the  Solution  Equation 

Before  attempting  a  numerical  solution  to  (10),  the  only  remaining  task  is  to 
specify  boundary  conditions  to  the  solution  area  detailed  in  Figure  6.  The  four 
boundaries  correspond  to  r  =  0,  r  =  /?wal],  emax,  and  u  =  0,  since  the  kinetic  energy  of  an 
electron  can  not  be  negative  as  the  total  energy  goes  to  zero. 


-V 


wall 


Figure  6  Diagram  of  typical  solution  area  including  the  four  boundaries. 
Boundaries  (2)-(4)  are  determined  at  run-time  by  the  input  conditions. 
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Two  of  the  boundaries  are  relatively  straightforward.  First,  for  the  cylindrical  systems 
usually  associated  with  the  positive  column,  the  radial  anisotropy,/,.,  must  be  zero  at 

r  =  0  in  order  to  remain  symmetric  around  the  axis. 


fr[r,£]  = 


1  d  7r  ? 


=  0 


(18) 


r=0 


{  K[r,£]  dr 

For  the  u  =  0  boundary  we  expand  (10)  and  take  the  limit  as  u  0  to  obtain  the  second 
condition.  This  leads  to 


(eoEzy  —  -e0Er[r]— 


/0[^£] 


=  0 


-e0V[r] 


(19) 


The  third  boundary  is  not  as  straightforward.  We  require  the  EEDF  to  go  to  zero 
as  £  — »  °o ,  and  this  is  in  fact  the  boundary  condition  chosen  by  Uhrlandt  and  Winkler  [7, 
pg  522].  Unfortunately,  the  numerical  solution  grid  must  end  at  some  finite  (preferably 
small)  upper  limit,  and  in  that  case  the  EEDF  is  not  zero.  We  ran  simulations  of  this 
boundary  condition  using  a  simple  model  that  had  an  analytic  answer.  These  simulations 
demonstrated  that  forcing  the  EEDF  to  zero  prematurely  introduced  significant  errors. 

The  condition  /0[r,£max  ]  =  0  is  not  valid  unless  the  model  is  run  to  extremely  high 

energies,  but  doing  so  exponentially  increases  the  computer  memory  requirements  and 
processing  time.  These  same  simulations  lead  to  the  resolution. 

The  analytic  simulations  showed  that  while  setting  the  last  energy  bin  to  zero 
introduced  significant  errors,  setting  it  to  a  fraction  of  the  previous  bin  did  not.  Choosing 
the  fraction  amounts  to  specifying  the  slope  of  the  EEDF  at  the  highest  energies  [15,  pg 

282],  and  allows  us  to  write  fQ[jmx  ]  =  Af0[jtmx  -1]  where  j  refers  to  the  energy  index 
and  A  the  fraction.  Early  versions  of  the  model  considered  only  elastic  collisions  with  a 
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constant  collisional  cross  section  and  no  radial  potential.  Under  these  assumptions,  (10) 
has  an  analytic  solution,  the  Druyvesteyn  distribution  (section  2.2).  This  analytic  result 
initially  provided  an  excellent  slope  for  the  top  boundary  condition.  However,  once  we 
included  inelastic  collisions  into  the  model,  the  distribution  became  distorted  enough 
from  a  Druyvesteyn  to  require  a  new  functional  form. 

In  order  to  successfully  utilize  the  fraction  A  as  an  upper  boundary  condition,  we 
needed  an  analytic  function  that  reflected  the  contribution  from  inelastic  collisions.  At 
the  highest  energies,  equation  (10)  takes  on  a  limiting  form  in  which  the  field  driven  flux 
toward  higher  energy  is  balanced  by  inelastic  collisional  losses.  This  equation  can  be 
solved  approximately,  providing  an  excellent  choice  for  A.  The  derivation  is  found  in 
Appendix  C.  With  it,  the  third  boundary  condition  becomes. 


MrJ™]  =  Af(l[r,jm„  -1] 

A  =  exp[-/3  (ujmm  -  M;max_! )] 


P  = 


K[r,emm]  H[r,em3 J 


(20) 


It  is  important  to  note  that  when  implemented  as  written  in  (20),  the  boundary 
condition  broke  down  near  the  wall.  As  r  approaches  /?wall,  the  kinetic  energy  u  at  emax 
decreases.  This  constriction  in  energy  space  violates  the  earlier  assumption  of  field 
driven  flux  balancing  collisional  losses.  To  overcome  the  problem,  we  used  umm[r=0] 

instead  of  Mmax[r]  when  calculating  both  p  and  A.  Another  possible  fix  would  be  to  find  a 

new  limiting  form  for  the  distribution  that  holds  true  near  the  wall.  But  this  would  also 
require  merging  the  two  boundary  conditions  at  an  intermediate  r.  As  it  stands,  the 
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current  method  works  well  for  the  parameters  of  interest,  and  is  superior  to  assuming 
7o[^emax]  =  0- 


The  last  boundary  condition  at  r  =  /?wall  is  also  difficult.  Electrons  with  sufficient 


total  energy  may  overcome  the  radial  potential  (e  >  -e0V[rwalt  ] )  and  be  “absorbed”  at  the 

wall  and  hence  lost  to  the  system.  Thus,  at  the  wall,  we  need  to  include  some  form  of 
loss  mechanism;  this  loss  will  also  balance  the  gains  from  ionization  to  produce  the 
steady  state.  Again,  initially  following  the  paper  of  Uhrlandt  and  Winkler  [7],  we  set  the 
radial  anisotropy  equal  to  a  proscribed  loss  function  at  the  wall. 

1  d  (fo Kan  > ^ exp[-au 2[rwall  ,£]] 


frKa,n£]  =  — 


(21) 


KKall>£]dr 

where  the  parameters  a  and  A  are  available  to  adjust  the  total  amount  of  loss  (see  section 
4.4). 


Uhrlandt  and  Winkler  experimented  with  a  number  of  different  wall  loss 
functions,  most  variations  of  the  exponential  form  in  (21)  [16].  They  came  to  the 
conclusion  that  except  for  the  region  very  near  to  the  wall,  the  solution  to  the  EEDF 
yields  essentially  the  same  results  independent  of  the  loss  function,  provided  that  the 
integrated  strength  is  the  same  [7,  pg  524], 

The  Uhrlandt  formalism  is  helpful  because  it  provides  inhomogeneous  terms  to 
resolve  our  solution,  lessening  the  reliance  on  Penning  ionization  (see  section  3.1.2). 
However  the  loss  function  form  is  only  loosely  physically  based.  The  search  for  a  more 
realistic  wall  loss  mechanism  led  to  recent  work  by  Alves,  et  al.  [17]  in  which  they 
presented  a  more  physically  based  boundary  condition  called  the  wall  loss  cone  (see 
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Figure  7).  Their  derivation  is,  in  turn,  an  extension  of  the  previous  work  by  both  Busch 
and  Kortshagen  [15]  and  Tsendin  and  Golubovskii  [18]. 


Uhrlandt  and  Winkler 
Loss  Function 

frKall’£]  =  gKall,£} 


-V, 


wall 


Alves,  et  al. 
Loss  Cone 


Figure  7  Comparison  of  the  two  types  of  wall  loss  boundary  conditions. 

Using  the  loss  cone  formalism,  electrons  must  have  enough  kinetic  energy 
perpendicular  to  the  wall  to  surmount  the  wall  potential.  An  additional  parameter,  the 
wall  loss  coefficient,  permits  only  a  fraction  of  the  penetrating  electrons  to  be  lost. 
Unlike  the  previous  wall  loss  development  by  Uhrlandt  and  Winkler,  this  formalism 
includes  a  region  known  as  the  plasma  sheath,  illustrated  in  Figure  7,  which  is  found  at 
the  tube  wall.  The  plasma  sheath  is  assumed  to  be  infinitely  thin,  and  its  main  property 
of  interest  is  the  drop  in  potential  that  occurs  across  it.  In  the  work  by  Alves,  the  sheath 
potential  is  taken  to  add  and  additional  -4  to  -8  eV  to  the  wall  potential  [17,  pg  899], 
We  wanted  to  consider  both  types  of  wall  loss  conditions  in  this  treatment  and 
therefore  attempted  to  integrate  both  into  the  code.  To  do  so,  we  chose  to  ignore  the 
addition  of  the  sheath  potential,  A (p ,  at  the  wall.  Since  our  treatment  can  accept  any 
arbitrary  radial  potential  (section  3.2.1),  we  are  always  able  to  add  an  effective  sheath 
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potential  by  modifying  the  input  near  the  wall.  This  allows  us  to  mimic  the  conditions 
considered  by  Alves,  while  including  the  Uhrlandt  formalism  as  well. 

Taking  A cp  =  0,  the  modified  loss  cone  boundary  condition  is. 


frKall’£]=- 


1 


—(f  \r  £l)=—  ^ _ —f\r  el 

Kir^te]drVoir'*'el)  2{\+Xr  ^ 


(22) 


where  %  is  the  reflection  coefficient,  which  runs  from  0  to  1.  Like  the  parameters  a  and 
A  in  (21),  x  adjusts  the  total  amount  of  loss  occurring  at  the  wall 

Unlike  the  Uhrlandt  formalism,  the  loss  cone  does  not  provide  an  inhomogeneous 
term  for  our  solution  set.  This  leaves  only  the  Penning  term  to  resolve  the  system.  But  in 
some  cases,  the  contribution  from  Penning  ionization  is  too  small  to  accurately  drive  the 
system,  resulting  in  an  unstable  solution.  When  this  happens  we  choose  to  accept  the 
homogeneous  equations,  and  resolve  the  system  by  assigning  the  EEDF  a  value  at  a 


particular  point  in  the  system,  e.g.  set  f0[rspecia,,£special]  =  1 .  This  allows  us  to  solve  the 


homogenous  system  rather  than  trying  to  force  a  suspect  inhomogeneous  solution. 

With  the  equations  developed,  and  the  boundary  conditions  formed,  we  are  ready 
to  find  a  numerical  solution  to  the  elliptic  differential  equation. 

3.2.  Numerical  Method  of  Solution 

The  numerical  method  used  to  solve  (10)  is  clearly  described  in  numerical 
analysis  textbooks  such  as  Burden  and  Faires  [19].  We  start  by  forming  a  grid  over  the 
solution  space.  This  grid  specifies  all  of  the  solution  points  within  the  space,  as  well  as 
the  boundary  points.  Next,  the  CBE  and  boundary  conditions  are  discretized.  This 
allows  us  to  form  an  equation  for  each  solution  point  involving  its  neighbors.  The 
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complete  set  of  equations  corresponds  to  a  linear  system  of  the  form  A  x  =  b ,  which  is 
solved  for  x  ( /0 )  using  Gaussian  Elimination. 

While  the  brief  outline  above  points  to  the  simplicity  of  the  numerical  method,  a 
number  of  challenges  arise  during  implementation  that  are  presented  in  the  following 
sections.  The  first  is  how  to  form  a  grid  over  the  solution  space  (shown  in  Figure  8).  The 
technique  used  to  do  this  is  covered  in  section  3.2.1.  In  forming  the  grid,  we  use  unequal 
spacing  in  both  the  radial  and  energy  directions.  This  unequal  spacing  requires  a  more 
complicated  set  of  differencing  equations  than  those  used  for  equal  grid  spacing.  Section 
3.2.2  details  this.  Finally,  the  linear  system  depends  on  a  matrix  A,  which  is  sparse 
banded.  While  simple  Gaussian  Elimination  can  be  used  to  solve  it,  a  more  sophisticated 
technique  discussed  in  section  3.2.3  is  employed  which  speeds  the  solution  and  reduces 
memory  requirements. 


Radius  (i) 

Figure  8  Diagram  of  basic  parameters  used  in  forming  the  solution  grid. 
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3.2.1.  Forming  the  Solution  Grid 


The  first  step  in  the  numerical  solution  is  to  create  a  grid  overlying  the  solution 
space.  However,  in  our  case,  the  shape  of  the  solution  space  depends  on  the  problem 
considered.  Figure  6  shows  that  the  lower  boundary  ( u  =  0 )  depends  on  the  radial 
potential  (-V[r]),  which  in  turn  depends  on  the  specific  parameters  of  the  discharge;  so 
the  gridding  technique  must  be  able  to  handle  a  “general”  input  potential. 

Although  using  a  uniform  mesh  is  desirable,  a  number  of  problems  immediately 
arise.  First,  the  grid  points  do  not  necessarily  intersect  the  u  =  0  boundary.  This  would 
force  us  to  interpolate  the  boundary  condition  and  solution  in  this  region.  Second,  the 
shape  of  the  potential  dictates  that  some  regions  require  a  finer  grid  than  others  do.  This 
is  especially  true  near  the  wall,  where  the  potential  can  change  rapidly.  These  difficulties 
convince  us  to  forgo  uniform  spacing  and  employ  a  non-equidistant  mesh. 

To  create  the  grid,  first  consider  Figure  8,  which  shows  the  relationship  between 
the  various  grid  parameters.  In  the  diagram,  NI  and  NJ  refer  to  the  maximum  bin  index 
in  the  radial  and  energy  directions  respectively.  SPEC  is  a  special  value  chosen  to 
optimize  the  density  of  grid  lines;  during  implementation  SPEC  is  normally  set  to  be 
NI/2.  Using  Figure  8  we  can  write  equations  which  relate  the  various  parameters. 


(23) 

(24) 

(25) 

(26) 
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At  this  point  the  equation  set  is  underdetermined.  In  order  to  resolve  the  system, 
we  choose  to  enforce  a  relationship  between  the  two  uniformly  spaced  regions.  We  force 
the  normalized  distance  between  0  and  SPEC  on  the  radial  axis  and  SPEC  and  NI  on  the 
energy  axis  to  be  the  same.  This  relationship  produces  equation  (27),  and  with  the  use  of 
(23)-(26)  allows  us  to  solve  for  Ae. 

SPEC  •  A r  _  (NI  -  SPEC)  •  Ac  (27) 

R  ~  -V 

,ywall  ywall 

-Vwall  +V[^rf~(NI  -  SPEC)  Ae] 

Ae  = - ^ -  (28) 

NI-SPEC 

With  these  six  equations  we  now  have  a  method  to  generate  the  grid: 

a)  Input  the  parameters  NI,  SPEC,  /?wall,  Vwal|,  emax.  Input  -V[r]  as  a  set  of  data 
points,  this  will  also  permit  calculation  of  the  inverse,  -V  ^r]. 

b)  Calculate  Ae  using  the  transcendental  equation  (28). 

e  +V 

c)  Find  Ar  using  (27)  and  NJ  =  - sss.  +  MI 

Ae 

d)  Create  part  of  the  radius  grid.  From  i=l  to  SPEC,  radius[i\=radius[i-l]+Ar 

e)  Create  the  energy  grid.  From  j=  1  to  SPEC,  energy\j]=  -V[radius\j]].  And  for 
j=SPEC+ 1  to  NJ,  energy\j] = energy\j- 1  ] + Ae. 

f)  Finish  the  radius  grid.  From  r=SPEC+ 1  to  NI,  radius[i]=  -VA[energy[i\\. 

After  completing  steps  a)  through  f),  the  grid  is  formed. 

As  mentioned  earlier,  each  intersection  on  the  grid  corresponds  to  either  a 
solution  point  or  boundary  point.  The  distinction  is  that  every  solution  point  has  a 
corresponding  equation  assigned  to  it,  while  the  boundary  points  are  used  only  to  form 
the  equations  of  their  neighbors.  In  this  treatment  we  decided  to  include  three  of  the 
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boundaries  in  the  solution  space,  reserving  the  top  boundary,  emax,  solely  as  a  boundary 
condition. 

Figure  9  shows  the  solution  grid  in  i-j  space;  notice  that  unlike  the  radius-energy 
grid  in  Figure  8,  the  spacing  is  uniform.  In  addition,  the  u  =  0  boundary  is  always  a 
straight  line,  regardless  of  the  input  potential  (-V[r]).  For  these  reasons,  it  is  often  easier 
to  think  of  the  solution  grid  in  i-j  space  rather  than  radius-energy. 


Figure  9  A  schematic  of  the  grid  in  i/j  space.  Solution  points,  m, 
are  numbered  sequentially  from  1  to  DIM. 

The  solution  points  are  numbered  sequentially  starting  at  the  origin,  m  =  1 ,  and 
proceeding  up  the  energy  axis  to  the  last  bin,  NJ- 1.  The  count  continues  in  the  next  radial 
bin,  starting  at  i  =  j,  and  continuing  to  the  last  point  when  i  =  NI  and  j  =  NJ- 1 .  This  last 
point  corresponds  to  the  total  number  of  equations  and  solution  points  in  the  grid,  and  is 
referred  to  as  DIM,  the  dimension  of  the  solution.  When  we  employ  the  finite 
differencing  scheme  (section  3.2.2)  we  will  need  to  refer  to  the  points  surrounding  every 
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solution  point.  In  the  energy  (/)  direction  we  can  simply  add  or  subtract  one  to  reference 
the  surrounding  solution  points.  In  the  radial  direction  however,  the  enumeration  is  not 
that  simple,  and  we  must  use  counting  variables  INCUP  and  INCDOWN.  There  are  also 
occasions  where  we  must  refer  to  points  two  radial  bins  away  from  a  solution  point;  in 
those  cases  we  use  the  variables  INC2UP  and  INC2DOWN.  With  the  help  of  Figure  9, 
we  derive  the  following  expressions  for  the  variables  in  the  counting  scheme, 


m  =  j-i  +  l  +  iNJ-1-^—^  (29) 

2 

2NJ-NI  +  2NINJ-NI2  ,,m 

DIM  = -  (jU) 

2 

INCUP  =  NJ  —  i  —  l  (31) 

INCDOWN  =  NJ-i  (32) 

INC2UP  =  2NJ  -  2i  -  3  (33) 

INC2DOWN  =  2NJ  -2i  +  l  (34) 


Now  that  each  solution  point  has  a  label,  m,  and  we  have  a  mechanism  to 
reference  the  surrounding  points,  we  are  ready  to  transform  the  equations  into  a  finite 
differencing  form. 

3.2.2.  Finite  Differencing  on  the  Non-Uniform  Grid 

In  order  to  write  an  equation  for  each  solution  point  in  the  grid,  equation  (10) 
must  be  transformed  from  an  analytic  PDE  into  a  differenced  equation.  In  the  case  of 
uniform  grid  spacing  this  means  using  well  known  formulae  such  as  the  three  point 
equation, 


dr 


- f[r,£ ] 


„  /Ki>£j-2/[/;,eJ]+/[r;_1,g>] 

h2 


(35) 
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which  approximates  the  second  order  derivative,  /[r;,e;],  on  a  uniform  grid  with 

step-size  h.  Since  the  grid  spacing  in  our  treatment  is  not  uniform,  we  require  slightly 
more  complicated  versions  of  (35). 

This  analysis  follows  the  example  of  Uhrlandt  and  Winkler  [7,  pg  528], 
employing  a  three  point  central-differencing  scheme  except  at  the  boundaries,  where  we 
invoke  either  forward  or  backward  differencing.  The  equations  are  derived  using  the 
Method  of  Undetermined  Coefficients,  outlined  in  most  numerical  analysis  textbooks. 
The  results  are  presented  below,  along  with  a  number  of  step  relations  used  to  simplify 
the  expressions. 


|-/[r,£]| 

dr 


U] 


d^ 

dr 


f[r,e] 


i-lj 


d  -  -|| 

Vrfl’-eX 


i+1  J 


-  h,.  =  B,  (a,.  •  (/w  j  "  Aj )  - 1 S(/M  j  -  Aj )) 

=  hlA.  =  Bt  (-a,  ■  (fMJ  -  ftJ )  -  (A  +  2)(fi_hj  -  4, )) 
^  hMJ  =  B,  ((2 + a,  )(/w j  -  /„, )  -  ~  fu )) 


(36) 

(37) 

(38) 


where  h  is  introduced  as  a  shorthand  notation,  and 

A  =V(rM-rw) 

ai=(ri-ri_l)/{rM-r)  (39) 

A  =1M 

0 

Similar  expressions  are  also  developed  for  — /[r,e]  where  5,  a,  and  /3  become  C,  ^ 

d£ 

and  S  respectively.  Equations  (36)-(38)  are  used  further  to  develop  expressions  for  the 
second  derivative.  In  the  case  of  central  differencing,  this  becomes 
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dr 


P[r,£] 


MlN 

dr 


£  B,  fe  (/>.,  Au  -  P^j )-  ft  (/UA-u  -  KA, ))  (40) 


M.j 


where  P[r,e ]  is  an  arbitrary  function.  Again,  very  similar  expressions  arise  for  the  second 
derivative  with  respect  to  energy. 

Using  (40)  and  the  corresponding  equation  for  energy,  we  can  transform  (10)  into 
finite  differenced  form.  At  the  grid  boundaries,  care  is  taken  to  use  either  a  forward  or 
backward  differencing  equation  such  as  (37)  or  (38),  and  to  substitute  in  the  appropriate 
boundary  condition.  The  complete  set  of  differenced  equations  is  rather  long,  so  it  is 
presented  in  Appendix  D. 

Once  the  CBE  is  transformed,  it  can  be  evaluated  at  each  solution  point.  Each 
solution  point  then  has  a  corresponding  equation,  and  we  are  ready  to  solve  the  resulting 
linear  system. 

3.2.3.  Solving  the  Sparse  Banded  Matrix 

Once  every  solution  point  in  the  grid  has  an  equation  associated  with  it,  the  entire 

problem  is  reduced  to  solving  a  linear  system  of  equations  of  the  form  A  x  =  b .  These 

systems  are  readily  solved  using  Gaussian  elimination,  and  many  algorithms  are  available 

which  implement  this  technique.  However,  simple  Gaussian  elimination  is  rather  slow 

since  the  number  of  operations  required  is  proportional  to  the  cube  of  the  number  of 
% 

matrix  elements.  Early  test  runs  on  a  Pentium  166  took  approximately  one  hour  to  solve 
a  system  with  only  50  radial  and  70  energy  (50  X  70)  bins.  Our  goal  was  to  do 
“production”  runs  with  approximately  300  X  500  bins,  which,  using  equation  (30), 
corresponds  to  an  A  matrix  with  over  1010  elements.  This  would  have  severely  taxed 
currently  available  super-computing  resources. 


38 


f 


a\ 
a01  a 


a 


*12 


*21 

31 


22 

*32 


a 


23 

ci'i'i  a 


42 


*33 
^43 

fl53  ~ 


*54 


a45 

a55  a 


a 


*65 

75 


(a) 


56 


^64  @6*  @66  @ 


86 


*67 
2  77 


■*87 


*78 


^21 

ya31 


@n 

@23 

@34 

@45 

a56 

@61 

@ 1 8 

@22 

@33 

@44 

@55 

U66 

@11 

@SS 

@32 

@'43 

@54 

@65 

ai6 

@81 

@42 

@53 

@64 

@15 

aS6 

(b) 

Figure  10  Sparse  Banded  Matrix  in  traditional  storage  (a)  and  reduced  storage  (b)  as 
required  by  DGBSV  routine.  Zero  elements  are  not  shown. 


But  the  A  matrix  is  special.  As  with  most  numerical  PDE  solutions,  the  A  matrix 
is  sparsely  populated  along  the  main  diagonal.  This  special  form  is  called  sparse  banded 
(see  Figure  10).  We  can  take  advantage  of  the  abundance  of  zeros  in  the  matrix  to  reduce 
both  the  program’s  memory  requirements  and  time  of  execution. 

The  obvious  way  to  decrease  the  memory  requirement  is  to  store  only  the  non¬ 
zero  bands  of  the  matrix.  This  requires  special  algorithms  to  create,  initialize,  and  access 
the  matrix  elements,  which  adds  time  to  the  overall  process,  but  the  reward  is  a  dramatic 
reduction  in  the  amount  of  memory  required.  In  the  case  of  a  300  X  500  grid,  the  number 
of  matrix  elements  held  in  storage  falls  from  1010  to  108  using  this  technique. 

Advance  knowledge  of  the  many  zeros  in  the  matrix  can  also  dramatically  reduce 
the  number  of  operations  needed  to  solve  the  system.  The  system  is  still  solved  by 
Gaussian  elimination,  except  the  algorithm  is  fine-tuned  to  speed  the  solution  as  much  as 
possible.  The  specific  algorithm  used  in  this  work  is  called  “DGBSV”  and  comes  from 
the  LAP  AC  set  of  routines  made  available  on  the  National  Institute  of  Standards  and 
Technology  Internet  site. 

DGBSV  requires  the  A  matrix  be  transformed  to  a  special  form  as  shown  in 
Figure  10.  Under  the  transformation,  each  diagonal  in  A  becomes  a  row  in  A'.  The  new 
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matrix  has  the  same  number  of  columns  as  A,  but  the  number  of  rows  is  reduced  to 
2NL+NU+1.  NL  is  the  number  of  lower  diagonals  and  NU  the  number  of  upper 
diagonals  in  A.  The  first  NL  rows  of  A'  are  filled  with  zeros  initially,  and  used  by  the 
routine  as  a  work  space. 

The  DGBSV  routine  comes  from  a  family  of  algorithms  created  to  solve  special 
matrices.  It  calls  two  other  functions,  DGBTRF  and  DGBTRS;  the  first  computes  the  LU 
decomposition  of  A',  while  the  second  performs  back  substitution  using  b  to  get  the 
answer.  This  presents  an  opportunity  to  solve  the  problem  quickly  with  multiple  forcing 
vectors.  Since  the  majority  of  the  execution  time  is  spent  finding  the  LU  decomposition 
of  A',  very  little  is  added  by  calling  DGBTRS  multiple  times.  It  also  opens  the  window 
to  efficiently  implement  a  time  dependent  solution  in  the  future. 

The  final  size  of  A',  and  hence  the  time  required  to  solve  it,  depends  directly  on 
the  number  of  diagonals  in  A,  and  the  number  of  diagonals  depends  directly  on  the 
differencing  scheme.  In  most  cases,  we  use  three  point  central  differencing,  which  results 
in  2NJ-X  diagonals  in  A.  However,  as  described  in  section  3.2.2,  the  lower  boundary 
requires  backward  differencing  in  the  radial  direction.  If  the  three  point  scheme  is  used, 
it  increases  the  number  of  diagonals  in  A  to  3NJ-3,  which  is  a  significant  increase  for 
large  NJ.  To  circumvent  this  problem,  we  chose  to  employ  only  two  point  differencing 
on  the  lower  boundary,  keeping  the  number  of  diagonals  in  A  at  2NJ-1.  While  the  two 
point  scheme  is  less  accurate  than  three  point,  it  is  only  used  on  the  lower  boundary, 
which  constitutes  a  small  fraction  of  the  solution  points,  and  testing  showed  that 
deviations  between  the  two  methods  appear  only  in  the  7th  significant  digit  of  the  result. 
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Using  the  DGBSV  algorithm  instead  of  straight  Gaussian  elimination 
dramatically  reduced  the  program’s  run  time.  The  50  X  70  bin  test  case,  which 
previously  took  about  an  hour  on  a  home  PC,  now  finishes  in  less  than  one  minute.  The 
largest  grids  (300  X  500)  are  run  at  the  Major  Shared  Resource  Center  on  a  SGI/Cray 
Origin  2000  computer.  In  those  cases,  the  program  requires  2  Gb  of  memory  and  takes 
about  30  minutes  to  finish  when  using  only  one  CPU.  The  execution  time  of  the  solution 
scales  approximately  as  the  square  of  the  number  of  solution  points. 
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IV.  Implementation  of  the  Strict  Kinetic  Solution 


In  this  chapter  we  step  through  four  distinct  phases  of  development  and  validation 
in  the  code.  At  each  step  we  gain  insight  into  the  electron  kinetics  within  the  discharge, 
as  well  as  establishing  the  capabilities  and  limitations  of  the  model.  The  first  phase  treats 
the  simple  case  of  a  homogeneous  system,  excluding  inelastic  collisions  and  assuming  a 
constant  elastic  collision  cross  section.  Under  these  restrictions,  we  have  an  analytic 
solution  to  validate  the  numerical  solution;  it  is  the  Druyvesteyn  distribution  discussed 
earlier.  Next,  we  add  a  non-trivial  radial  potential  to  the  model,  comparing  the  numerical 
results  to  an  analogous  analytic  solution.  In  the  third  phase,  we  begin  to  add  source  and 
loss  terms.  While  not  physically  representative  collision  terms,  these  point  sources  and 
losses  demonstrate  some  of  the  capabilities  of  the  model  and  elucidate  behavior  typical  of 
the  local  and  nonlocal  regimes.  Finally,  we  complete  the  development,  adding  the 
realistic  collision  terms  as  well  as  wall  loss  to  the  system. 

4.1.  Phase  I:  Radially  Homogeneous  System 

This  first  version  of  the  model  neglects  inelastic  collisions  and  the  influence  of  a 
radial  electric  field,  and  assumes  a  constant  elastic  cross  section  (£>[u]=<2).  With  these 
simplifications,  (10)  becomes 


ru 


r  dr  3 K[r,e]dr 


d  d  u(e  E  )2  d  ,  d  A 

/«,+di(i^ir87/«)+^(G[r’£l/»)=0 


(41) 


where  u  equals  e,  since  the  radial  potential  is  assumed  zero.  As  there  are  no  volumetric 
sources  in  the  system  and  a  steady  state  solution  is  sought,  there  can  be  no  loses,  and  thus 
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the  wall  boundary  condition  is  set  to  perfect  reflection.  In  other  words  the  flux  of 


electrons  at  the  wall  is  set  to  zero, 


l 


r=Rwall 


=  0. 


Because  the  system  is  homogeneous  by  design,  we  know  that  the  local 
approximation  is  satisfied  (see  section  2.2).  As  noted  previously,  assuming  a  constant 
elastic  cross  section  and  no  inelastic  collisions  the  CBE  has  an  analytic  solution;  it  is  the 
Druyvesteyn  distribution  given  in  equation  (3)  and  shown  in  Figure  3b.  Recall  that  in  the 

local  approximation,  the  radial  term  (liLf — Hi — 7  \)  is  assumed  negligible  and 

rdr  3K[r,£]dr  0 

dropped,  whereas  in  the  numerical  model  it  is  included  but,  if  implemented  correctly,  still 
yields  the  local  result. 

We  ran  the  model  on  a  uniform  grid  with  NI  =  10  and  NJ  =  100 ,  using  only  ten 
radial  bins  since  the  system  was  radially  homogeneous.  In  order  to  compare  the  model 
output  (/0[r,e])  to  the  analytic  result  (F[r,e]),  we  normalized  f0  at  each  solution  point, 

i.e.  F[r,e]  =  f0[r,e]/ne[r] .  Comparing  the  numerical  model  output  to  the  analytic  result 

at  each  solution  point,  we  obtained  relative  errors  between  the  two  no  greater  than  1%. 
The  next  step  was  to  add  a  spatially  varying  radial  potential. 

4.2.  Phase  II:  Include  Radial  Potential 

Using  the  gridding  technique  presented  in  section  3.2.1,  we  added  a  general  radial 
potential  to  the  model.  At  this  point,  we  still  neglected  all  inelastic  collisions;  thus 
equation  (41)  remains  valid.  However,  with  the  addition  of  the  potential,  the  system  is  no 
longer  radially  homogeneous,  and  the  local  approximation  does  not  necessarily  hold. 
Although  the  system  is  no  longer  guaranteed  to  satisfy  the  local  approximation,  we  can 
force  it  to  this  regime  with  a  judicious  choice  of  input  parameters.  Since  we  had  already 
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validated  Phase  I  against  the  local  approximation,  we  chose  here  to  compare  to  the 
nonlocal  approximation. 

As  noted  in  section  2.3,  we  can,  under  certain  assumptions,  analytically  solve  (41) 
using  the  nonlocal  approximation.  With  no  inelastic  collisions,  a  constant  elastic  cross 
section,  and  a  quadratic  form  for  the  radial  potential,  we  arrive  at  the  modified 
Druyvesteyn  distribution,  equation  (7).  It  is  important  to  note  that  this  1-D  derivation 
assumes  that  all  electrons  have  a  total  energy  less  than  <p[/?wall],  which  is  effectively  the 
same  as  assuming  an  infinite  wall  potential.  Because  our  treatment  is  2-D,  it  obviously 
includes  a  finite  limit  on  (p  at  r  =  /?wall.  To  accommodate  this  difference,  we  reduced  emax 

to  9>[/?waU]. 

In  the  previous  phase  (section  4.1)  the  system  was  homogeneous,  so  we  knew  it 
satisfied  the  local  approximation;  but  in  this  case  we  had  to  ensure  the  system  was  in  the 
nonlocal  regime.  This  meant  setting  the  PR  value  such  that  the  electron  energy  relaxation 
length  was  much  greater  than  the  tube  radius,  the  opposite  requirement  of  equation  (4). 

We  ran  the  model  for  hypothetical  gas  at  a  pressure  of  0.1  Torr  and  r  =  1cm 
( PR  =  0.1 ).  As  required,  this  PR  is  much  less  than  the  3  Torr-cm  established  on  the 
right-hand-side  of  equation  (5).  The  addition  of  the  radial  potential  forced  us  to  employ 
the  non-uniform  grid  of  section  3.2.1,  and  we  chose  to  set  NI  =  100  and  NJ  =  105  in  an 
effort  to  approximate  the  energy  mesh  density  of  Phase  I.  Using  these  parameters,  we 
again  obtained  relative  errors  between  the  analytic  and  numerical  results  of  less  than  1% 
throughout  the  volume. 

Having  added  the  radial  potential  into  the  system,  we  moved  on  to  include 
volumetric  sources  and  losses. 
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4.3.  Phase  III:  Add  Source  and  Loss  Terms 


With  the  addition  of  inelastic  collisions,  the  model  would  be  complete.  But 
before  addressing  these  complicated  terms,  we  investigated  how  the  model  reacts  to  very 
simple  losses  and  gains.  As  in  the  previous  section,  we  continued  to  neglect  inelastic 
collisions,  assume  a  constant  elastic  collision  cross  section,  and  use  a  quadratic  radial 
potential.  Equation  (41)  is  used  to  describe  the  system,  along  with  a  “total  reflection” 
wall  boundary  condition. 

We  started  by  adding  a  simple  point  source  to  the  volume  at  a  given  radius  and 
energy.  In  this  case  the  source  term  is  an  inhomogeneous  addition  to  (41),  e.g.  + 1014 ,  at 
one  particular  solution  point,  r',  ef.  Since  there  are  no  losses  at  the  boundaries,  a 
volumetric  loss  is  included  to  maintain  a  steady  state.  As  with  the  source  term,  the  loss  is 
established  at  a  single  computational  point,  in  this  example  r,  £".  While  the  source  term 
is  inhomogeneous,  the  loss  is  proportional  to  the  local  population.  Global  particle 
conservation  then  fixes  the  population  at  the  loss  point  such  that  1014  =  Af0[r',£"] , 

where  A  is  the  proscribed  factor  of  proportionality. 

Based  on  physical  reasoning  alone,  we  expect  that  at  the  source  and  loss  points 
the  EEDF  will  display  “deviations”  from  the  unperturbed  Druyvesteyn  form  (Figure  3b 
page  9).  Since  electrons  are  being  added  at  the  source  point  (r,  ef),  the  distribution 
function  there  should  be  slightly  higher  than  neighboring  radial  bins.  Conversely,  at  the 
loss  point  (r,  £"),  where  electrons  are  removed  from  the  system,  the  distribution  should 
show  a  deficit.  We  find  in  fact  that  this  behavior  is  highly  dependent  on  the  pressure 
regime  of  the  system. 
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Running  the  model  at  different  pressures,  we  noticed  a  manifestation  of  local 
versus  nonlocal  behavior.  In  the  local  (high-pressure)  case,  the  electron  energy  relaxation 
length  is  short;  hence  disturbances  in  the  distribution  function  are  confined  to  the  radial 
location  of  the  source/loss.  In  the  nonlocal  regime  the  opposite  is  true,  and  any 
disturbance  in  the  distribution  is  propagated  to  all  radii.  Figure  1 1  illustrates  these 
differences. 


Figure  1 1  Log/linear  plots  of  the  EEDF  versus  total  energy  at  five  neighboring  radial  bins 
demonstrating  the  effect  of  a  volumetric  point  gain/loss  in  the  local  (a)  and  nonlocal  (b) 
regimes.  The  plots  are  vertically  offset  in  order  to  distinguish  differences  between  them. 

In  the  figure  above  we  see  the  results  of  the  volumetric  point  gain/loss,  and  how 
its  effect  is  propagated  in  radius.  In  both  cases  the  solid  line  identifies  the  distribution  at 
/,  the  radius  at  which  the  loss  occurs,  while  on  either  side  of  it,  the  dashed  lines  show  the 
distributions  of  four  neighboring  radial  bins.  The  curves  are  vertically  offset  in  order  to 
distinguish  the  differences  between  them.  For  the  local  case  (a),  the  distribution  is 
reduced  dramatically  at  the  point  of  the  loss,  but  only  two  radial  bins  away  has  recovered 
almost  completed.  In  contrast,  the  nonlocal  regime  shows  less  depletion  of  the 
distribution  at  r  and  the  deficit  is  apparent  at  all  r.  The  log  scale  of  the  figure  masks  the 
effect  of  the  point  source  in  both  cases. 
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The  results  discussed  above  are  different  than  those  obtained  in  a  one-dimensional 


solution,  as  shown  in  Figure  12.  In  a  1-D  model,  such  as  the  local  or  nonlocal 
approximation,  the  distribution  function  is  unable  to  recover  from  a  point  loss,  and  takes 
on  a  reduced  form  at  all  energies  greater  than  e".  In  the  examples  of  Figure  11,  both 
distributions  recover  at  higher  energies  because  electrons  are  able  to  flow  around  the  loss 
point.  These  electrons  feed  the  high-energy  tail  of  the  distribution,  maintaining  its 
overall  shape.  In  a  1-D  model,  there  is  no  way  for  the  electrons  to  flow  around  the  loss 
point;  thus  the  distribution  never  recovers  from  the  loss. 

If  point  losses  are  included  at  all  radial  bins,  generating  a  line-of-losses  in  the  e-r 
space,  we  obtain  results  analogous  to  the  1-D  models.  As  explained  above,  the  reason  the 
distribution  functions  in  Figure  1 1  recover  at  higher  energies  is  because  electrons  are  able 
to  flow  around  the  loss  point,  feeding  the  tail.  Extending  our  model  from  a  point  loss  to  a 
line-of-losses  impedes  this  flow.  Once  the  loss  extends  across  all  radii,  there  is  no  longer 
a  path  to  the  high-energy  tail,  and  the  distribution  takes  on  a  reduced  form  identical  to  a 
1-D  result. 
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Legend 


Figure  12  Log/linear  plot  of  the  EEDF  versus  energy  for  the  case  of  a  point  source  and 
an  uninterrupted  line-of-losses  within  the  2-D  strict  model.  Model  output  is  also 
equivalent  to  the  result  of  a  1-D  model  with  a  single  point  loss.  Dashed  lines  show 
undisturbed  Druyvesteyn  distributions. 

Figure  12  displays  model  output  for  the  case  of  a  line-of-losses.  Although  only 
one  radial  location  is  shown,  the  effect  is  the  same  at  all  radii.  In  this  figure,  the  solid 
line  indicates  the  model  output,  while  the  dashed  lines  are  undisturbed  Druyvesteyn 
distributions.  The  lower  curve  has  been  shifted  down  to  align  it  with  the  high-energy  tail 
of  the  output.  Notice  that  below  25  eV  the  output  is  unaffected  by  the  loss-line.  Where 
the  loss  occurs,  there  is  a  steep  drop,  then  at  energies  above  30  eV  the  distribution 
recovers  to  the  Druyvesteyn  form,  although  reduced  in  magnitude.  This  is  exactly  the 
same  behavior  seen  for  a  point  loss  in  a  1-D  model. 

These  limiting  examples  demonstrate  the  ability  of  the  model  to  react  to  losses 
and  gains  within  the  system.  They  also  serve  once  again  to  illustrate  differences  between 
the  local  and  nonlocal  regimes. 
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4.4.  Phase  IV:  Completion  of  the  Model 


After  gaining  insight  using  artificial  sources  and  losses,  we  were  ready  to  add 
realistic  collision  terms  and  complete  the  model.  This  section  discusses  implementation 
details  of  the  final  solution. 

By  including  inelastic  collisions,  the  model  moves  from  equation  (41)  to  the  full 
CBE,  equation  (10),  and  the  associated  difference  equations  in  Appendix  D.  The  addition 
of  volumetric  ionization  requires  us  to  choose  one  of  the  wall  loss  formalisms,  (21)  or 
(22),  in  order  to  assure  a  steady  state  solution.  Having  added  this  final  functionality,  the 
model  itself  is  complete.  However,  in  order  to  accurately  model  physical  systems  we 
must  also  specify  physically  representative  input  parameters. 

A  self-consistent  model  is  one  that  uses  only  externally  determined  parameters 
such  as  the  pressure  in  order  to  reach  a  solution.  In  the  case  of  the  positive  column,  a 
completely  self-consistent  model  would  require  only:  constituent  gas  data,  gas  partial 
pressures,  axial  field  strength,  and  axial  current.  But  the  model,  in  its  present  form,  is  not 
self-consistent  regarding  three  internal  variables:  the  radial  potential,  excited  state 
number  density,  and  total  wall  loss.  Therefore,  in  addition  to  the  externally  determined 
parameters,  these  three  internal  variables  must  be  specified  as  well. 

In  the  present  model,  values  for  the  radial  potential  are  either  assumed  or 
measured  experimentally.  In  this  thesis  we  use  assumed  forms  exclusively,  usually 
quadratic  or  cubic.  Experimental  data,  if  available,  are  preferred,  but  their  use  is  limited 
to  the  specific  parameters  of  the  experiment.  Another  option  is  to  first  run  a  self- 
consistent  nonlocal  model,  and  then  use  the  calculated  potential  as  input  to  the  strict 
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solution.  While  this  method  is  only  an  approximation,  in  most  cases  it  is  probably  more 
accurate  than  simply  assuming  a  form. 

A  similar  situation  exists  with  the  excited  state  number  densities,  which  are  used 
in  the  calculation  of  de-excitation  collisions  and  Penning  ionization  (see  section  3.1.2). 
These  collisions  can  have  a  significant  impact  on  the  EEDF,  if  the  excited  state  density  is 
a  large  fraction  of  the  total  density.  In  most  physically  realistic  cases  however,  we  have 
seen  that  excited  state  neutrals  represent  a  small  fraction  of  the  overall  neutral  population. 
We  therefore  feel  confident  using  educated  estimates  for  these  values,  rather  than  pursue 
a  self-consistent  treatment.  Future  extension  of  the  solution  to  a  quasi-self-consistent 
model  is  discussed  in  chapter  VII. 

Having  estimated  values  for  the  radial  potential  and  excited  number  density,  we 
are  left  with  one  remaining  internal  variable  -  the  total  wall  loss.  While  the  wall 
boundary  conditions  (21)  and  (22)  proscribe  the  form  of  the  loss,  both  have  parameters  to 
adjust  the  total  integrated  loss.  We  use  this  parameter  to  adjust  the  total  number  of 
electrons  in  the  system,  i.e.  adjust  the  magnitude  of  the  EEDF.  Increasing  the  total  loss 
forces  the  number  of  electrons  to  increase  as  well,  in  order  to  maintain  the  steady  state.  If 
measurements  of  the  axial  discharge  current  are  available,  this  adjustment  allows  us  to 
match  the  computed  axial  current,  equation  (53),  to  the  experimental  value.  Without 
these  experimental  measurements  we  are  forced  to  estimate  the  value  of  the  total  wall 
loss  and  adjust  the  wall  loss  parameter  accordingly. 
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V.  Model  Validation 


To  validate  the  completed  model,  we  compare  its  output  to  previously  published 
work  as  well  as  to  available  approximate  models.  Since  our  solution  formalism  closely 
follows  the  development  of  Uhrlandt  and  Winkler,  comparison  to  their  work  should  yield 
nearly  identical  results.  In  the  limit  of  high  and  low  pressure,  we  can  verify  our  results 
against  those  of  the  local  and  nonlocal  approximations,  using  code  previously  developed 
by  Bennett  [6],  Finally,  we  validate  the  model  against  a  totally  different  method,  namely 
the  Monte  Carlo  simulation  technique. 

5.1.  Exposition  of  Model  Output  and  Comparison  to  Uhrlandt  and  Winkler  Results 

Output  from  our  strict  solution  is  compared  to  the  work  of  Uhrlandt  and  Winkler 
[7]  for  a  neon  discharge.  The  parameters  for  the  discharge  are  a  column  radius  of 
^waii  =1-7  cm,  a  gas  pressure  of  0.62  Torr,  which  corresponds  to  a  density  of 

N0  =  2.0xl016cm~3 ,  and  an  axial  current  of  Iz  =  10  mA.  The  axial  field  strength  is  2.1 
V/cm  and  the  radial  potential  has  a  minimum  at  the  wall  of  V[r]  =  -21  V.  We  chose  to 
use  the  Uhrlandt  wall  loss  condition  (equation  (21)),  setting  a  =  -0.1  and  A  =  1.05xl012 
in  order  to  satisfy  the  proscribed  axial  current.  For  inelastic  collisions  we  considered  two 
excited  states  of  neon,  a  metastable  ( m )  and  resonance  (r)  level,  each  with  population 
densities  of  2x10"  cm-3 .  Further  data  for  those  two  states  are  shown  in  Table  1. 
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Table  1  Electron/Atom  and  Atom/Atom  Collision  Data  Used  in  Neon  Discharge 
Calculations.  Notation:  m  -  metastable,  r -  resonance,  o  -  ground  state. 


Electron/Atom  Collision  Process 

Energy  loss/gain  (eV) 

Excitation 

Ne0+e  — >Nem+e 

-16.62 

Ne0  +  e~  — >  Ner  +  e~ 

-16.67 

De-excitation 

Nem  +e~  — >  Ne0  +e~ 

16.62 

Ner  +e~  — >  Ne0  +e~ 

16.67 

Ground  State 
Ionization 

Ne0  +  e~  Ne+  +  2e~ 

-21.56 

Atom/ Atom  Collision  Process 

Rate  Coefficient  (cm3  s1) 

Penning  Ionization 

Nem  +  Nem  — >  Ne0  +  Ne+  +  e~ 

Zn  3.4xlO”10 

Ner  +  Ner  — >  Ne0  +  Ne+  +  e~ 

z22  3.4  xlO-10 

Nem  +  Ner  — >  Ne0  +  Ne+  +  e~ 

Zj2  6.8xlO-10 

It  is  important  to  note  a  number  of  differences  between  the  two  implementations 
of  this  strict  solution.  Both  solutions  admit  a  generalized  potential,  but  Uhrlandt  used  an 
experimentally  derived  form,  while  we  chose  a  simple  cubic  since  we  didn’t  have  access 
to  the  same  data.  In  addition,  Uhrlandt  performed  a  side  calculation  to  estimate  the 
excited  state  number  densities.  Although  our  treatment  does  not  include  this  side 
calculation,  we  based  our  excited  densities  on  the  published  results  of  Uhrlandt.  The 
input  densities  do  not  match  exactly,  because  our  model  assumes  they  are  uniformly 
distributed,  while  Uhrlandt  allows  them  to  vary  with  r.  Even  so  the  results  are 
remarkably  similar. 

Figure  13  displays  the  three  scalar  components  of  the  EEDF,  f0[r,u ],  /r[r,w],  and 
fz[r,u],  at  two  different  radii,  r  =  0.6  and  r  =  1.2  cm.  While  examples  of  the  isotropic 
distribution, /0,  have  been  shown  previously,  this  is  the  first  example  of  the  radial  and 
axial  components.  Note  that  at  all  radii,  the  axial  component  of  the  distribution  is 
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negative,  reflecting  the  action  of  the  axial  electric  field,  directed  in  +  z ,  to  bias  the 
electron’s  velocity  in  the  -  z  direction.  The  radial  distribution,  on  the  other  hand, 
undergoes  a  sign  change.  At  lower  energies,  a  positive  ft  corresponds  to  a  radial  flux 
toward  the  wall,  while  at  higher  energies,/,,  turns  negative  for  much  of  the  solution  area, 
indicating  flow  away  from  the  wall.  At  some  point  this  high  energy  flow  must  reverse 
sign  in  order  to  satisfy  the  wall  loss  boundary  condition;  this  is  discussed  further  in 
Chapter  VI. 


Figure  13  Log/linear  plot  of  components  of  the  EEDF,  for  the  neon  discharge  at  two  radii, 
calculated  using  the  strict  solution:  isotropic  -f0[r,u],  radial  -  ft[r,u],  and  axial  -fz[r,u]. 

A  comparison  of  Figure  13  with  that  found  in  the  Uhrlandt  paper  [7,  pg  537] 
shows  excellent  agreement  for  all  three  components.  At  both  radial  locations,  the  peak 
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magnitude  of/0  matches  exactly  with  those  in  the  Uhrlandt  paper,  and  the  overall  forms 
of  the  distributions  are  identical.  The  two  anisotropic  components  are  also  very  similar; 
in  Figure  13,  fr  changes  sign  at  energies  of  10.5  and  15  eV,  while  in  Uhrlandt’s  paper  the 
reversal  occurs  at  1 1  and  15  eV.  The  anisotropic  components  also  show  many  of  the 
same  small-scale  variations  in  both  figures. 

While  a  graphical  comparison  of  our  EEDF  to  Uhrlandt’s  work  is  satisfying,  more 
rigorous  tests  exist.  From  the  distribution  function,  we  derive  a  host  of  macroscopic 
observables,  which  are  often  the  true  quantities  of  interest  (see  Appendix  A).  Space  does 
not  permit  us  to  display  all  of  the  results,  but  we  obtained  excellent  agreement  with 
Uhrlandt  for  quantities  such  as  the  electron  number  density,  mean  energy,  and  integrated 
currents.  However,  the  most  rigorous  test  of  the  model  comes  from  looking  at  particle 
and  energy  balances.  As  presented  in  Appendix  A,  the  particle  balance,  equation  (54), 
equates  the  particles  lost  at  a  given  radius  (divergence  of  particle  current)  to  the  particles 
gained  through  ionization  sources  (ground  and  Penning  ionization  rate).  Similarly,  the 
energy  balance,  equation  (55),  compares  the  net  energy  leaving  a  given  radius 
(divergence  of  the  energy  current)  to  that  gained  and  lost  from  collision  terms  as  well  as 
the  action  of  the  axial  and  radial  electric  fields. 

In  Figure  14,  we  see  a  plot  of  the  particle  balance  for  the  neon  case  considered.  If 
the  model  achieves  complete  particle  balance,  the  sum  of  the  two  ionization  rates  (ground 
state  ionization  and  Penning  ionization)  will  exactly  equal  the  divergence  of  the  particle 
current  at  each  radial  location.  As  seen  in  the  figure,  the  total  ionization  rate  for  neon  is 
everywhere  slightly  less  than  the  net  divergence.  However,  we  obtained  relative  errors  in 
the  particle  balance  less  than  2%  for  the  majority  of  the  tube,  which  is  the  same  as  that 
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reported  by  Uhrlandt.  In  the  region  between  0.9  and  1.2  cm,  we  see  a  greater  deviation  of 
approximately  12%  in  the  balance  that  is  as  yet  unexplained. 


Figure  14  Terms  in  the  electron  particle  balance  as  a  function  of 
radius  for  the  case  of  a  neon  discharge. 

If  we  compare  Figure  14  to  the  equivalent  figure  in  Uhrlandt’ s  paper  [7,  pg  542], 
we  notice  immediately  the  difference  in  the  Penning  ionization  treatment.  As  stated 
earlier,  Uhrlandt  included  a  radial  dependence  in  the  excited  state  species  concentration, 
with  a  maximum  on  axis  that  decays  to  zero  at  the  wall.  Because  of  this,  his  Penning 
term  reflects  the  same  radial  dependence,  contributing  zero  at  the  wall.  In  our  model 
however,  the  excited  state  species  are  uniformly  distributed,  which  results  in  an  equal 
Penning  contribution  at  all  radii. 
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Figure  15  Important  terms  in  the  local  energy  balance  as  a  function  of  radius  for 
the  neon  discharge.  Not  shown  are  elastic  collisional  losses  and  inelastic 
collisional  gains,  which  are  insignificant  in  this  case. 


The  radial  energy  balance  is  shown  in  Figure  15.  As  in  the  case  of  the  particle 
balance,  the  divergence  of  the  energy  current  should  equal  the  sum  of  the  various  energy 
gains  and  losses  at  each  radial  location.  It  is  important  to  note  that  the  signs  of  the 
divergence  and  inelastic  loss  terms  are  reversed  from  equation  (55)  in  an  effort  to 
condense  the  figure.  As  plotted,  the  divergence,  axial  heating,  and  radial  cooling  sum  to 
give  the  inelastic  loss.  In  addition,  the  magnitude  of  two  terms,  the  losses  due  to  elastic 
collisions  (Lelastic)  and  the  energy  gains  from  collisions  (Gc),  are  too  small  to  appear  in  the 
figure. 

Figure  15  also  corresponds  to  a  similar  figure  in  Uhrlandt’s  paper  [7,  pg  543].  In 
his  results,  Uhrlandt  quotes  relative  errors  in  the  energy  balance  of  less  than  0.1%.  For 
the  neon  case  shown,  this  treatment  achieved  errors  less  than  0.5%  for  the  majority  of  the 


56 


tube  radius.  Near  r  =  1.2  cm  the  divergence  approaches  zero,  and  the  relative  error 
increases  dramatically;  however,  this  is  not  considered  an  indication  of  any  error  in  the 
model. 

The  figure  above  demonstrates  why  energy  transport  is  observed  from  the  outer  to 
inner  parts  of  the  neon  discharge.  Near  the  axis,  inelastic  losses  outweigh  axial  heating, 
but  the  terms  reverse  dominance  at  approximately  r  =  0.75  cm.  Excess  energy  added 
near  the  wall  thus  moves  toward  the  axis  to  fill  the  deficit  there.  Also  interesting  is  the 
effect  of  the  radial  field.  A  net  radial  current  flows  outward  to  supply  the  wall  loss;  this 
flow  is  along  the  radial  field,  which  (for  negatively  charged  electrons)  leads  to  a  small 
radial  cooling. 

5.2.  Verification  Against  the  Local  and  Nonlocal  Approximations 

Local  and  nonlocal  approximate  models  present  another  opportunity  to  validate 
our  solution.  While  numerous  published  results  exist,  we  have  the  capability  to  use  “in- 
house”  versions  of  these  approximations  developed  by  Bennett  [6].  In  this  way  we  can 
ensure  the  input  parameters  are  exactly  the  same.  In  addition,  we  are  able  to  compare  the 
model  output  directly  rather  than  examine  diagrams  from  a  journal  article. 

In  these  tests  we  want  to  validate  the  strict  model,  not  the  approximations.  With 
this  in  mind,  the  input  parameters  are  chosen  to  create  the  conditions  for  which  the 
respective  approximations  are  most  valid.  As  an  example,  neither  approximation 
includes  a  mechanism  for  wall  loss;  thus  for  this  validation,  we  employ  the  wall  loss  cone 
formalism,  with  the  reflection  coefficient  set  to  1  -  total  reflection.  Without  wall  loss,  we 
are  also  forced  to  neglect  ionization  sources,  in  order  to  satisfy  the  steady  state 
requirement. 
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The  model  validation  was  made  using  neon  as  the  neutral  gas,  with  the  same  input 
parameters  as  used  in  the  Uhrlandt  comparison  above,  unless  specified  otherwise.  We 
chose  input  values  to  guarantee  that  the  system  fell  into  the  appropriate  regime,  local  or 
nonlocal,  when  comparing  to  the  respective  models.  Based  on  the  results  of  equation  (5), 
we  consider  the  neon  column  local  at  a  PR  =  50  and  nonlocal  at  PR  =  0. 1 .  Keeping  the 
ratio  of  E/N  constant,  this  leads  to  the  following  input  parameters:  local  comparison  - 
P  =  29.4  Torr,  Ez  =  99.6  V/cm;  nonlocal  comparison  -  P  =  0.062  Torr,  Ez  =  0.21  V/cm. 

One  difficulty  that  exists  in  the  comparisons  is  that  the  local  and  nonlocal  models 
use  uniform  gridding,  while  the  strict  model  uses  a  non-uniform  grid.  Thus  the  output 
data  do  not  coincide,  and  we  can  not  calculate  the  relative  error  between  specific  solution 
points  in  the  models.  To  overcome  this,  the  output  from  both  models  is  interpolated 
using  cubic  splines.  The  interpolation  polynomials  are  then  sampled  uniformly,  using  the 
resulting  data  set  to  calculate  the  relative  error. 

Another  consequence  of  the  difference  in  gridding  is  how  the  cross  section  data  is 
sampled.  The  elastic  collision  cross  section  of  neon  has  a  feature  at  low  energies  (the 
Ramsauer  minimum)  that  influences  the  shape  of  the  distribution  function.  The  non- 
uniform  grid  technique  (section  3.2.1)  results  in  high  resolution  of  this  low  energy 
feature,  while  for  the  same  number  of  bins,  the  uniform  grid  could  not  resolve  it.  To 
correct  this  deficiency,  the  local  and  nonlocal  models  were  run  on  a  mesh  with  twice  the 
number  of  energy  bins  as  the  strict  solution. 

Using  the  input  parameters  and  methodology  discussed  above,  the  models  were 
compared  over  a  range  of  30  eV  and  at  three  radial  locations  (r  =  0,  0.5  and  1  cm).  The 
results  from  the  nonlocal  case  were  uniformly  impressive.  The  models  deviated  by  no 


58 


more  than  4%  in  relative  error,  while  the  distribution  itself  spanned  over  4  orders  of 
magnitude.  In  the  local  comparison,  the  results  were  equally  impressive  at  the  first  two 
radial  locations,  with  the  relative  error  again  below  4%.  But  at  r  =  1  cm,  the  wall 
boundary,  the  models  deviated  significantly.  We  believe  this  is  due  to  the  rapid  decrease 
in  the  radial  potential  near  the  wall,  which  causes  a  breakdown  in  the  local  assumption  of 
no  radial  variations  over  the  electron  energy  relaxation  length  (section  2.2).  As  such, 
these  errors  reflect  on  the  local  approximation,  and  not  the  implementation  of  the  strict 
solution. 

5.3.  Comparison  to  Monte  Carlo  Methods  Using  Argon 

In  the  previous  two  sections  we  compared  our  strict  kinetic  solution  results  to 
those  of  other  kinetic  solutions.  As  a  final  check  on  the  model,  we  look  to  a  completely 
different  technique  for  verification,  Monte  Carlo  particle  simulation.  Specifically  the 
strict  solution  is  run  under  the  conditions  reported  by  Kortshagen  et  al.  [9],  in  which  they 
investigate  the  EEDF  of  an  argon  discharge. 

As  with  the  Uhrlandt  comparison  in  section  5.1,  we  again  try  to  match  the  input 
parameters  to  those  stated  in  the  paper.  For  the  radial  potential,  Kortshagen  assumed  a 
quadratic  form.  At  the  wall  boundary,  he  employed  a  Monte  Carlo  version  of  the  loss 
cone  boundary  condition,  assuming  A (p  =  8  eV  and  total  absorption  at  the  wall,  %  =  0.  As 
explained  in  section  3.1.3.,  our  treatment  of  the  loss  cone  formalism  differs  slightly, 
incorporating  the  potential  jump,  A (p ,  as  a  steep  increase  in  the  radial  potential  just  prior 
to  the  wall,  rather  than  a  step  function  at  r  =  /?wall.  Finally,  Kortshagen  made  no  mention 
of  including  the  effects  of  Penning  ionization  or  collisions  of  the  2nd  kind,  thus  these 
features  were  turned  off  for  the  model  comparison. 
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Figure  16  Comparison  (a)  of  Monte  Carlo  calculation  [9]  (left)  and  strict 
solution,  present  work,  (right)  for  an  argon  column  at  three  different  radii. 

Plotted  is  the  EEDF  as  a  function  of  total  energy. 

We  compared  the  two  models  for  two  cases,  (a)  and  (b);  the  relevant  input 
parameters  are  given  in  Table  2.  Figure  16  presents  the  results  for  the  first  comparison: 
on  the  left  a  photocopy  from  the  Kortshagen  paper,  on  the  right  the  output  from  our  strict 
model.  A  similar  presentation  is  made  in  Figure  17  for  the  second  set  of  input 
parameters. 


Table  2  Input  parameters  for  model  comparisons  between  Monte  Carlo 
and  strict  solutions  in  an  argon  column. 


Comparison  (a) 

Comparison  (b) 

Number  density 

N  =  lxl016cm“3 

/V  =  3xl016cnT3 

Pressure 

P  =  0.31  Torr 

P  =  0.93  Torr 

Axial  Field 

Ez  =  6.53  V/cm 

Ez  =  12.4  V/cm 

Sheath  Potential 

Vsheath  =  -8.97  V 

Sheath  = -10-31V 

Wall  Potential 

Vwall  = -16.97  V 

Vwall  = -18.31V 

60 


Figure  17  Comparison  (b)  of  Monte  Carlo  calculation  [9]  (left)  and  strict 
solution,  present  work,  (right)  for  an  argon  column  at  three  different  radii. 

Plotted  is  the  EEDF  as  a  function  of  total  energy. 

Although  the  two  models  use  vastly  different  approaches  to  simulate  the  argon 
column,  the  results  are  remarkably  similar.  In  Figure  16,  the  plots  at  all  three  radii  follow 
almost  exactly  those  of  Kortshagen.  In  the  second  run,  Figure  17,  a  slight  deviation  is 
noted  for  the  EEDF  at  r  =  0.96  cm,  but  this  is  attributable  to  our  incorporation  of  the 
sheath  potential  increase,  A (p ,  into  the  radial  potential  rather  than  as  a  step  increase  at  the 
wall. 


These  results,  combined  with  those  of  the  previous  two  sections,  demonstrate  the 
accuracy  of  this  strict  kinetic  model  as  compared  to  other  well-established  solutions. 
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VI.  Investigation  of  Special  Topics 


As  mentioned  initially  in  the  introduction,  the  main  purpose  of  this  thesis  was  to 
develop  and  validate  the  strict  solution  method.  Continuing  the  theme,  this  chapter  looks 
more  closely  at  some  of  the  characteristics  of  the  strict  solution.  In  doing  so,  we  further 
our  understanding  of  the  solution  and  at  the  same  time  reveal  interesting  properties  of  the 
electrons  kinetics  in  a  glow  discharge. 

6.1.  Current  Flow  within  the  Solution  Space 

An  exciting  consequence  of  the  strict  solution  is  that  we  can  model  radial  current 
flow  within  the  plasma.  Since  the  approximate  methods  are  one-dimensional  and 
disregard  wall  loss,  they  result  in  no  net  radial  flow.  In  this  section,  we  revisit  the  results 
of  the  neon  discharge  discussed  previously  in  section  5.1,  using  the  same  input 
parameters. 

As  listed  in  Appendix  A,  the  radial  and  axial  particle  currents  are  derived  from  the 
anisotropic  components  of  the  EEDF  by  integrating  the  particle  flux  over  energy 
(equation  (48)).  Considering  the  differential  flux,  we  can  map  electron  flow  within  the 
solution  grid.  Figure  18  presents  a  vector  field  plot  of  the  electron  flux,  or  differential 
particle  current,  within  the  neon  solution  area.  The  figure  is  constructed  by  combining 
the  radial  and  axial  flux  components  into  a  vector;  all  of  the  vectors  are  then  scaled  to  the 
same  length.  Unfortunately,  this  scaling  removes  all  information  about  the  magnitude  of 
the  vectors.  Note  that,  although  it  is  not  apparent  from  the  figure,  the  differential  flux 
becomes  negligibly  small  in  the  upper  left  comer  of  the  solution  grid. 
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Figure  18  Vector  diagram  of  electron  flux  within  the  neon  discharge  solution  area. 

Three  mechanisms  contribute  to  create  a  particle  flux  circuit  in  the  solution  space 
shown  in  Figure  18:  particle  diffusion,  the  action  of  the  radial  electric  field,  and  inelastic 
collisions.  Starting  near  the  origin  of  the  diagram,  we  see  that  diffusion  dominates, 
driving  low  energy  electrons  toward  the  wall,  along  (but  against  the  influence  of)  the 
radial  electric  field.  Diffusion  therefore  produces  a  small  cooling  of  the  electrons.  At  the 
edge  of  the  tube,  some  high-energy  electrons  exit  the  system  via  the  wall  boundary 
condition.  The  rest  reverse  course  and  head  back  toward  the  axis,  driven  now  by  the 
radial  electric  field.  To  complete  the  circuit,  inelastic  collisions  scatter  high-energy 
electrons  back  to  lower  energies.  The  majority  of  this  scattering  takes  place  near  the  axis, 
because,  for  a  given  total  energy,  electrons  there  have  the  greatest  kinetic  energy.  Of 
course,  against  the  backdrop  of  all  this  radial  flow  is  the  constant  action  of  the  axial  field, 
driving  electrons  to  higher  energies. 
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Figure  19  Integrated  current  flows  within  the  neon  discharge, 
plotted  as  a  function  of  the  tube  radius. 

A  complementary  view.  Figure  19,  displays  the  currents  for  the  same  neon 
column;  these  currents  are  found  by  integrating  the  differential  flux  shown  in  Figure  18 
over  energy  (equations  (48)  -  (52)).  The  net  radial  current,  jfamcl\  is  the  sum  of  two 
components:  the  flow  due  to  diffusion,  j^l‘n ,  and  that  due  to  the  radial  electric  field, 

jfieulcle  •  As  described  earlier,  the  diffusion  current  is  comprised  of  low  energy  electrons 

moving  toward  the  wall,  while  high-energy  electrons  make  up  the  field  driven  current. 
Notice  that  the  diffusion  and  field  driven  currents  are  nearly  equal  and  opposite,  although 
when  summed,  they  result  in  the  net  positive  particle  current  ( j?ar“cle),  which  supports  the 
wall  loss.  While  the  net  current  is  positive,  most  of  the  high-energy  electrons  are  moving 
towards  the  axis,  resulting  in  a  negative  energy  current  ( jermrgy ).  Only  near  the  wall, 
where  all  of  the  current  becomes  positive,  does  the  energy  current  change  sign. 
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It  is  important  to  reiterate  that  only  through  the  2-D  strict  solution  are  we  able  to 
investigate  these  radial  currents. 

6.2.  Applicability  of  the  Approximate  Methods 

The  strict  kinetic  method  provides  a  robust  and  flexible  model  of  the  EEDF,  but  it 
is  computationally  demanding.  Because  of  this,  we  might  be  tempted  to  employ  one  of 
the  approximate  kinetic  methods,  if  they  are  valid.  As  discussed  in  Chapter  n,  the  local 
(nonlocal)  approximation  is  applicable  when  the  electron  energy  relaxation  length  is 
much  less  (greater)  than  the  tube  dimensions.  While  expressions  such  as  (5)  attempt  to 
formalize  this  requirement,  they  do  not  provide  absolute  guidance  concerning  when  the 
approximate  solutions  may  be  used.  Ingold  suggests  a  general  range  of  applicability, 
which  is  listed  in  Table  3. 


Table  3  Range  of  validity  for  models  of  interest  [8,  pg  5943], 


Method 

PR<  1 

1<PR<10 

10</W<100 

Strict  Kinetic 

Yes 

Yes 

Yes 

Nonlocal  Approximation 

Yes 

? 

No 

Local  Approximation 

No 

No 

? 

We  know  that  neither  approximate  solution  will  be  accurate  for  all  r,  because 

neither  includes  a  wall  loss.  The  nonlocal  approximation  assumes  that  =  0 

dr 

everywhere,  which  is  violated  for  any  nonzero  wall  loss  since  fr  °c  -^2-  >  Furthermore,  in 

dr 

the  local  case,  the  radial  potential  varies  rapidly  in  the  vicinity  of  the  wall  (Figure  2,  page 
8),  which  violates  the  assumption  that  the  electron  relaxes  in  energy  space  before 
conditions  change  radially.  The  question  therefore,  should  not  be  “when  are  the 
approximate  methods  valid?”,  but  “how  accurate  can  we  expect  them  to  be?” 
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We  have  already  shown  in  section  5.2  that  in  their  regimes  of  applicability,  both 
the  local  and  nonlocal  approximations  yield  results  within  4%  of  the  strict  solution  over 
at  least  half  the  tube  radius.  However,  in  that  validation  phase  the  rather  unrealistic 
condition  of  zero  wall  loss  was  chosen  to  minimize  the  differences  in  the  models;  we  now 
wish  to  revisit  the  question  using  a  more  stringent  evaluation. 

Because  our  strict  solution  accurately  spans  both  high  and  low  pressure  regimes, 
we  can  use  it  alone  to  determine  the  validity  of  the  nonlocal  and  local  approximations. 

For  instance,  in  the  nonlocal  approximation,  the  EEDF  is  a  function  of  total  energy  only. 
Thus,  to  test  the  “nonlocality”  of  the  system,  we  simply  compare  the  values  of  /0  at 

different  radii,  keeping  e  constant.  If  the  system  is  completely  nonlocal,  then  the  relative 
difference  between  f0[rv  ej  and  f0[r2,  ej  will  be  zero.  Deviations  from  nonlocality 

correspond  to  errors  in  the  approximation  and  hence,  errors  inherent  in  the  nonlocal 
model.  We  can  therefore  use  the  deviations  to  gauge  the  accuracy  of  the  approximation. 

As  in  section  5.2,  we  ran  the  tests  for  a  neon  column.  The  following  parameters 
were  used:  /?wall  =  1  cm,  a  quadratic  potential  with  <p[/?walI]  =  20  V,  and  a  loss  cone 

condition  with  %  =  0.5.  The  resulting  data  are  presented  in  Table  4.  It  shows  the 
maximum  relative  difference  between /0  on  axis  and  at  three  other  radial  values;  this 
deviation  is  also  mapped  as  a  function  of  PR  value. 

Table  4  Deviation  from  nonlocality  (relative  percent)  for  a  neon  column 
as  a  function  of  PR  value  and  radius. 


PR  (Torr-cm) 

r  =  0.5  i?wall 

r  =  0.75  i?wall 

r  ^wall 

0.01 

<  1% 

<  1% 

i% 

0.05 

1% 

3% 

7% 

0.1 

3% 

6% 

14% 

0.5 

2% 

10% 

46% 
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The  results  in  Table  4  show  that  the  applicability  of  the  nonlocal  approximation 
depends  not  only  on  PR,  but  radius  as  well.  Note  the  disparity  between  these  results  and 
the  general  guidance  in  Table  3.  According  to  Ingold,  the  nonlocal  approximation  is 
valid  for  PR  =  0.5,  but  in  our  tests  at  r  =  /?wa„  the  relative  error  is  46%.  In  all  cases  the 
maximum  deviations  from  nonlocality  occur  in  the  high  energy  tail  of  the  distribution 
function,  reflecting  in  part  the  effect  of  the  wall  loss  boundary  condition. 

It  is  also  interesting  to  compare  how  other  gases  react  to  the  test  of  nonlocality. 
For  neon,  equation  (5)  gives  an  approximate  PR  value  of  17  Torr-cm.  Argon,  on  the 
other  hand,  results  in  a  value  of  ~  6  Torr-cm.  We  therefore  expect  that  for  the  same  PR 
value,  an  argon  column  will  have  greater  deviations  from  nonlocality  than  neon.  Table  5 
shows  the  results  of  the  argon  column  using  the  same  parameters  as  for  neon.  As 
expected,  the  deviations  are  much  greater  in  each  case. 


Table  5  Deviation  from  nonlocality  (relative  percent)  for  an  argon 
column  as  a  function  of  PR  value  and  radius. 


PR  (Torr-cm) 

r  =  0.5  i?walI 

r  =  0.75  /fwa]1 

r  ^wall 

0.01 

1% 

3% 

7% 

0.05 

4% 

8% 

30% 

0.1 

14% 

22% 

35% 

0.5 

96% 

99% 

100% 

We  can  make  similar  comparisons  in  the  local  regime.  If  the  column  satisfies  the 
local  approximation  exactly,  the  normalized  EEDF  will  be  the  same  at  all  r  regardless  of 
the  radial  potential,  i.e.  f0[rl,u]/ne[rl]  =  f0[r2,u]/ne[r2].  Unfortunately,  in  this  case  we 
are  comparing/0  at  constant  u,  which  means  that  the  grid  points  do  not  coincide.  To 
work  around  this  problem,  we  employ  the  same  technique  as  in  section  5.2,  interpolating 
the  data  and  then  re-sampling  it  uniformly  to  calculate  the  relative  errors.  This  technique 
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worked  well  for  the  neon  column,  but  failed  for  argon.  The  failure  stemmed  from  the 
inability  of  the  spline  interpolation  routine  to  adequately  fit  the  model  output  over  the  40- 
decade  drop  in  the  distribution  function.  The  spline  routine  returned  a  highly  oscillatory 
fit  that  resulted  in  extremely  large  relative  errors.  So  we  ran  the  comparison  for  neon 
only,  using  the  same  input  parameters  as  above. 


Table  6  Deviation  from  locality  (relative  percent)  for  a  neon 
column  as  a  function  of  PR  value  and  radius. 


PR  (Torr-cm) 

r  =  0.5  flwall 

r  =  0.75  i?wall 

r  =  Kali 

10 

19% 

55% 

350% 

50 

7% 

23% 

145% 

The  most  important  observation  from  the  data  above  is  how  poor  the  local 
approximation  is  at  r  =  /?wall.  This  again  relates  to  the  fact  that  near  the  wall,  two 
assumptions  are  violated.  First,  as  mentioned  in  section  2.2,  at  the  wall  the  radial  electric 
field  is  not  negligible  compared  to  the  axial  field.  And  second,  the  potential  varies 
rapidly  in  this  area,  violating  the  assumption  of  radial  homogeneity  over  the  scale  of  an 
electron  mean  free  path. 

After  including  the  effects  of  wall  loss,  it  is  apparent  from  the  data  above  that  the 
approximate  methods  are  limited  to  a  range  of  pressure,  as  well  as  a  fraction  of  the  tube 
radius. 

6.3.  Effect  of  the  Wall  Loss  Boundary  Condition 

As  first  explained  in  section  3.1.3,  the  wall  loss  boundary  condition  is  an 
extremely  important,  yet  almost  arbitrary  restriction  on  the  solution.  We  choose  simple 
expressions  to  represent  the  loss,  masking  the  true  nature  of  this  problem,  which  involves 
complicated  physical  processes  including  ion/electron  recombination,  absorption,  and 
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reflection  [20].  While  the  wall  loss  cone  has  some  physical  basis,  both  it  and  the  wall 
loss  function  largely  neglect  these  complicated  processes. 

In  his  paper,  Uhrlandt  mentions  that  he  tested  a  number  of  different  forms  for  the 
loss  function,  including  various  exponential  and  trigonometric  functions,  but  found  that 
in  general  they  “yielded  almost  the  same  results  for  the  EYDF  except  for  the  region  very 
near  to  the  wall”  [7].  In  our  studies,  we  found  that  while  this  statement  was  true  for  the 
two  cases  Uhrlandt  presented,  it  is  not  true  in  general. 


Figure  20  Comparison  of  the  on-axis  EEDF  for  the  wall  loss  cone  (dashed)  and  wall  loss 
function  (solid)  at  two  PR  values:  (a)  PR  =  1.0  and  (b)  PR  =  0.1. 

Both  of  Uhrlandt’ s  studies  considered  gases  in  which  the  PR  value  was  greater 
than  one.  As  shown  in  Figure  20(a),  for  PR  ~  1 ,  the  two  boundary  conditions  are  indeed 
very  similar  far  away  from  the  wall.  However,  as  we  should  expect  from  the  earlier 
discussions  of  the  nonlocal  regime,  at  low  PR,  the  EEDF  is  a  function  of  total  energy 
only,  and  any  differences  that  exist  at  the  wall  are  reflected  on  axis  as  well.  Part  (b)  of 
the  figure  shows  that  at  PR  ~  0.1  the  solutions  diverge  significantly,  even  on  axis.  We 
thus  need  the  “insulation”  provided  by  high  PR  to  confine  our  uncertainties  in  boundary 
condition  to  the  region  near  the  wall. 
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Within  Uhrlandt’s  wall  loss  function,  it  is  interesting  to  look  more  closely  at  how 
different  choices  for  the  parameter  a  (equation  (21),  page  29)  effect  the  form  of  the 
EEDF.  As  explained  earlier,  we  used  a  =  -0. 1  throughout  this  thesis,  as  Uhrlandt  did  in 
his  treatment,  but  we  never  gave  a  physical  justification  for  this  choice.  The  parameter  A, 
on  the  other  hand,  is  chosen  at  run-time  so  the  output  matches  experimental  data  (section 
4.4). 


Figure  21  Isotropic  component  of  the  EEDF  at  r  =  RwaU  for  three 
different  values  of  the  wall  loss  function  parameter  a. 

Figure  21  illustrates  how  three  different  values  of  the  parameter  a  effect  the  shape 
of  the  isotropic  component  of  the  EEDF,  for  PR  =  1 .  As  the  magnitude  of  a  increases, 
the  exponential  loss  function  drops  more  quickly  to  zero  with  increasing  energy;  thus  to 
support  a  fixed  loss,  the  majority  of  the  loss  occurs  at  low  energies.  On  the  other  hand,  as 
a  decreases  the  loss  function  diminishes  less  rapidly  with  energy,  and  hence  more  of  the 
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loss  occurs  at  high  energies.  This  effect  is  reflected  in  the  sharp  dip  seen  in  the  EEDF  at 
low  energies  for  a  =  -1.0 ,  as  opposed  to  the  sustained  decrease  seen  for  a  =  -0.05 . 

Since  the  wall  loss  function  is  arbitrary,  and  only  loosely  physically  based,  we 
really  have  no  basis  for  one  choice  of  a  over  another.  In  the  end,  it  comes  down  to 
aesthetics.  From  the  figure,  we  see  that  when  a  =  -0.1 ,  the  model  produces  an  EEDF 
which  is  smooth  and  relatively  flat  on  a  log  scale.  In  addition,  this  particular  choice 
produces  an  EEDF  similar  to  that  given  by  the  approximate  methods;  it  thus  has  the 
shape  we  have  come  to  expect.  Other  choices  for  a,  especially  a  =  -1 .0 ,  produce 
distortions,  which  without  any  physical  justification  to  use  them,  we  choose  not  to  accept. 
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VII.  Conclusions  and  Recommendations 


7.1.  Conclusions 

The  primary  objective  of  this  thesis  was  to  implement  and  test  a  strict  kinetic 
solution  to  the  CBE  of  a  glow  discharge.  We  have  achieved  this  goal,  and  thus  fashioned 
a  powerful  tool  to  explore  the  electron  kinetics  of  weakly  ionized  non-equilibrium  gases. 
In  the  process,  we  developed  numerical  techniques  that  are  directly  applicable  to  the 
solution  of  related  PDEs  and  gained  insight  into  the  applicability  of  other  approximate 
kinetic  models.  Furthermore,  we  treated  important  mechanisms  such  as  radial  particle 
and  energy  transport,  which  are  neglected  in  the  approximations. 

The  general  numerical  techniques  discussed  in  chapter  III  comprise  a  robust  set  of 
tools  applicable  to  a  variety  of  numerical  problems.  In  section  3.2.1,  we  presented  a 
method  to  generate  non-uniform  grids  over  a  variable  solution  area.  This  gridding 
scheme  required  us  to  develop  non-uniform  differencing  equations,  as  well  as  special 
counting  variables  to  access  the  points  in  the  solution  grid.  In  section  3.2.3,  we  discussed 
our  use  of  publicly  available  algorithms  to  solve  sparse  banded  matrices  and  detailed  an 
efficient  memory  storage  scheme  for  the  model.  Our  treatment  of  the  resulting  system  of 
equations  reduced  the  computational  load  significantly  over  “brute-force”  methods. 

During  the  course  of  the  investigation  we  frequently  validated  our  model  against 
approximate  kinetic  methods.  We  saw  that,  given  the  proper  assumptions,  the  strict 
solution  matched  very  closely  both  analytic  (section  4.1  and  4.2)  and  numerical  results 
(section  5.2)  of  the  local  and  nonlocal  approximations.  In  later  investigations  however, 
we  recognized  that  for  realistic  systems,  the  approximate  methods  failed  to  adequately 
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describe  the  entire  solution  area,  due  in  large  part  to  their  lack  of  a  wall  loss  condition 
(section  6.2).  We  concluded  therefore,  that  only  a  strict  solution  provides  a  complete 
description  of  the  positive  column,  especially  at  the  operating  pressures  of  today’s 
physical  devices. 

Throughout  the  validation  and  investigation  phases  of  the  thesis,  we  repeatedly 
demonstrated  the  necessity  of  employing  a  two-dimensional  solution.  It  was  only 
through  this  strict  solution  that  we  were  able  to  properly  account  for  radial  current  flows 
and  energy  transport,  and  explain  such  phenomenon  as  radial  electron  cooling  and  the 
negative  energy  current  in  a  neon  column. 

After  devoting  much  of  this  thesis  to  the  development  and  validation  of  the  strict 
kinetic  model,  we  have  only  begun  to  demonstrate  its  capabilities.  The  material 
presented  here  accounts  for  a  mere  fraction  of  the  possible  permutations.  In  fact,  we 
considered  only  single  species  atomic  gases  throughout  the  treatment,  but  our  solution 
also  accommodates  gas  mixtures  and  molecular  gases6.  In  its  present  form,  practical 
applications  for  the  solution  are  enormous.  The  model  is  ready  to  be  used  to  optimize  the 
kinetic  design  of  lighting,  electric  discharge  lasers,  and  plasma  processing  devices.  But 
development  should  not  end  here.  As  mentioned  previously  in  the  Introduction,  this 
research  marks  a  first  generation  effort  toward  the  future  exploitation  of  plasmas  in 
modifying  aircraft  flight  characteristics.  With  this  goal  in  mind,  a  number  of  model 
extensions  are  recommended  as  future  research  projects. 


6  We  are  able  to  treat  molecular  gases  by  including  the  vibrational  and  rotational  energy  level  cross- 
sections  and  increasing  the  energy  grid  density  in  order  to  resolve  the  finer  details  of  these  interactions. 
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7.2.  Recommendations  for  Future  Work 

The  most  obvious  extension  of  the  model  is  to  make  it  time  dependent,  rather  than 
steady  state.  This  modification  only  involves  the  inclusion  of  one  additional  term  into  the 
CBE.  And  as  noted  in  section  3.2.3,  our  current  numerical  technique  is  well  suited  to 
employ  this  time  dependence.  By  adding  time  dependence,  the  solution  could  then  be 
used  to  investigate  wave  propagation  as  well  as  plasma  instabilities  in  the  positive 
column. 

A  more  radical  variation  would  be  to  switch  from  a  study  of  radial  to  axial 
kinetics.  The  present  solution  models  the  axially  homogenous  positive  column  region  of 
the  glow  discharge,  but  a  change  in  coordinates  would  allow  investigation  of  axial 
inhomogeneities.  One  practical  use  for  this  is  the  study  of  shocks  in  a  discharge  tube. 
Shocks  can  be  modeled  by  rapid  variations  in  the  constituent  number  density  and 
resulting  space  charge  field.  Our  current  solution  accepts  a  general  potential,  which  in 
the  new  problem  would  allow  for  a  jump  in  the  electric  field  strength;  however,  the 
model  does  not  include  a  mechanism  to  vary  the  number  density  spatially.  We  could  add 
a  spatial  dependence  to  the  density  in  two  ways,  either  as  a  constantly  varying  parameter 
like  the  potential,  or  as  a  step  function,  taking  on  two  different  values,  one  on  either  side 
of  the  shock.  The  most  difficult  part  of  the  transition  would  be  to  specify  appropriate 
conditions  at  the  two  spatial  boundaries.  Unlike  the  cylindrical  case,  there  are  no 
symmetries  to  take  advantage  of  at  these  boundaries.  One  possible  technique  would  be  to 
specify  the  absolute  axial  flux  at  both  boundaries,  similar  to  Uhrlandt’s  wall  loss 
function. 
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Further  modifications  to  the  code  would  allow  us  to  investigate  an  even  wider 
class  of  physical  problems.  As  an  example,  the  solution  method  could  be  applied  to 
studies  in  the  space  environment.  As  with  shocks  in  a  discharge  tube,  after  a  simple 
change  of  coordinates  and  boundary  conditions,  we  should  be  able  to  employ  the  strict 
solution  technique  to  solve  for  the  EEDF  and  resulting  macroscopic  quantities.  Of  course 
in  this  case,  we  would  also  need  to  add  a  magnetic  field  to  the  acceleration  terms. 

During  the  development  of  our  solution,  we  noted  a  number  of  limitations  to  the 
model,  some  of  which  could  also  be  the  subject  of  future  research.  The  two  most 
important  involve  the  excited  state  number  densities  and  the  radial  potential.  In  both 
cases,  the  current  solution  method  does  not  solve  for  these  parameters  self-consistently; 
rather  it  leaves  them  as  variables  to  be  determined  through  experiment  or  supplementary 
analysis. 

Uhrlandt  and  Winkler  have  employed  a  quasi-self-consistent  technique  that  solves 
for  the  excited  state  number  densities.  Their  technique  uses  a  complementary  calculation 
to  the  strict  solution,  adjusting  the  excited  state  densities  iteratively  after  a  rate  balance 
analysis  [7].  In  order  to  add  this  functionality  to  our  model,  the  first  step  would  be  to 
make  the  excited  state  number  density  a  function  of  radius.  With  the  radial  variation 
included,  the  next  step  would  be  to  add  Uhrlandt’ s  rate  balance  analysis  as  a  new 
component  to  the  existing  code.  Once  included,  the  modified  code  would  require  an 
iterative  approach,  starting  with  an  assumed  form  for  the  number  densities  and  updating 
the  fraction  of  excited  neutrals  after  each  iteration,  until  the  computed  axial  discharge 
current  converges  on  an  experimental  or  expected  value. 
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The  radial  potential  also  presents  an  opportunity  for  future  work.  To  date  no  one 
has  performed  a  self-consistent  strict  solution  for  the  radial  potential,  although  Bennett 
[6]  has  done  so  for  the  nonlocal  approximation.  In  order  to  incorporate  this  calculation, 
we  would  have  to  simultaneously  determine  a  radial  distribution  for  the  ions  as  well  as 
the  electrons.  The  resulting  distributions  would  then  be  used  to  solve  Poisson’s  equation 
for  the  radial  potential.  While  we  could  theoretically  write  a  CBE  for  the  ions  and  solve 
it  directly,  Bennett  has  shown  that  for  the  conditions  of  interest  a  far  more  simple 
technique  is  to  employ  a  fluid  analysis  for  the  ions,  while  solving  for  the  electron 
distribution  directly.  This  is  analogous  to  the  method  Uhrlandt  used  in  his  treatment  of 
excited  state  number  densities.  After  including  the  fluid  calculation  into  the  code,  the 
solution  proceeds  iteratively,  starting  with  an  assumed  form  for  the  radial  potential  and 
updating  it  after  each  calculation,  until  it  converges  on  a  self-similar  form. 

In  both  these  self-consistent  treatments,  as  well  as  any  time  dependent  approach,  a 
limiting  factor  on  the  model’s  usefulness  is  program  run-time.  We  mentioned  in  section 
3.2.3,  that  the  current  code  actually  calls  two  algorithms.  The  first,  and  most 
computationally  intensive,  finds  the  LU  decomposition  of  the  input  matrix  (A),  while  the 
second  uses  that  result  along  with  the  forcing  vector  to  solve  the  system  of  equations.  If, 
between  iterations  or  time-steps,  the  A  matrix  does  not  change  much,  we  can  bypass  the 
decomposition  and  solve  the  system  quickly.  In  practice  we  might  update  the  matrix  only 
once  every  5  or  10  iterations  or  time-steps,  vastly  reducing  the  overall  run-time.  Even  so, 
both  the  self-consistent  and  time  dependant  problems  will  require  much  more  computing 
power  than  used  in  this  thesis. 
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Appendix  A  Macroscopic  Quantities  Derived  from  the  Distribution  Function 


While  the  result  of  solving  (10)  is  the  isotropic  distribution  function,  as  a  practical 
matter  we  are  more  concerned  with  macroscopic  quantities  such  as  the  axial  current, 
ionization  rate,  and  electron  mean  energy.  Not  only  are  these  the  quantities  used  to 
design  functioning  equipment,  but  they  constitute  the  easiest  values  to  compare  our 
results  to  other  published  works. 

All  of  the  macroscopic  quantities  are  derived  from  an  energy  space  integration 
over  one  of  the  three  distribution  functions  f0,  fr,  or  fz,  where  the  anisotropic  distributions 
are  given  by  equations  (15)  and  (16)  respectively.  Averaging  over  /0  leads  to  the 
electron  density,  ne[r],  the  mean  energy,  ue[r],  and  the  mean  electron  collision  rates 

Pu  M: 


ne[r]  =  jf0[r,u]uU2du 


(42) 


ue[r]  =  j  f0[r,u]uV2du  Jne[r ] 


PkreSS[r]  =  J— J  Nk  QresSWfo[r,u]u  du 

V  me  0 


(43) 

(44) 


The  mean  ionization  collision  rate  is  found  by  substituting  Qi°niza,ion  jnt0  (44),  with 

analogous  substitutions  to  determine  the  excitation  and  de-excitation  rates. 

Energy  loss  densities  associated  with  elastic  and  inelastic  processes  are  given  by, 


Lelastic[r]  =  X^J—  f  NkQeklas[u]f0[r,u]  u2du 


k  Mk  \l  m 


e  0 


j  inelastic  r  -i  _  V  V  excitation  p  excitation  .  V*  ionization  p  ionization  r  n 

^  iri-  2j2juu  +2juk  n  in 


(45) 

(46) 


k  l 
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and  associated  energy  gains  determined  from 


Gc[r ]  =  ^  [r]  +  ££  up™ning  Pf?nmg  [r]  (47) 


k  l 


k  k* 


Using  the  anisotropic  components  of  the  distribution,  we  obtain  particle  and 
energy  current  densities 


7r“W  =  ;j— \f,[r,u]udu 
0 


1  2 


"  ‘  J  3  l)m. 


jr‘’lr]  =  U—]fxlr,u]u2du 

3  \me  0 

where  x  refers  to  either  r  or  z.  The  radial  particle  current  is  composed  of  two  parts, 

■  particle  r  i  _  *  parricte  r  l  ,  •  partite  r  -i 

Jr  L'  J  J  diffusion  L '  J  '  J  field  driven  LA  J 


(48) 

(49) 


(50) 

(51) 

(52) 


where  AT[r,e]  is  given  by  (83).  The  total  axial  discharge  current  is  given  by 

Rwall 

lz=-2ne0  jjz[r]rdr 

o 

Finally,  the  solution  is  checked  for  self-consistency  using  particle  and  energy 
balances.  The  local  particle  balance  equation  for  electrons  is 


(53) 


~(r  jp;rMe  [r])=  X  P“m  M  +  XX  Pkpk*nnins  [r]  (54) 


k  k* 


And  the  local  electron  energy  balance  is  given  by 


(r  je;ergy[r])=  -e0Er[r]jrlde  ~ e0Ez [r]jpzar,ide  +  Gc[r ] - LMc[r] - L'"[r]  (55) 


78 


Appendix  B  Derivation  of  Strict  Solution  Equations 


This  appendix  presents  a  more  thorough  derivation  of  the  equations  given  in 
section  3.1.1.  It  is  based  entirely  on  the  previous  developments  of  Holstein  [10,  12]  and 
Allis  [13,  21].  As  in  section  3.1,  the  derivation  begins  with  the  collisional  Boltzmann 
equation, 

1 (56) 

where  /  =  f[r,v,t]  is  the  electron  velocity  distribution  function  (EVDF);  the  other  terms 
are  described  more  fully  on  page  17. 


The  treatment  in  this  thesis  considers  a  steady  state  solution,  in  which  case 


6  — 

and  axial  component,  thus  a  =  — -E . 

m 

Under  the  conditions  described  in  section  3.1.1,  we  make  the  assumption  that  the 
EVDF  is  nearly  isotropic,  allowing  an  expansion  in  spherical  harmonics.  Keeping  only 
the  first  two  terms  gives 


fir,  v]  =  XP,(c°s[0])/((r,  v) 
/ 


V 


(57) 


where  /0  is  the  isotropic  part  of  the  distribution,  and  /,[r,v]  =  frer  +  fzez  is  the 
anisotropic  part. 
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The  next  step  is  to  substitute  the  results  of  (57)  into  (56)  and  simplify  the  result. 


Considering  each  piece  separately,  we  have  for  the  radial  flux  term 

v-Vr/  =  |v|Vr(cos[0]/)=|v|Vr^cos[0]/;P;[cos[0]] 

i 

Using  the  identity, 

cos[fl]P,[cos[fl]]  =  (/  +  1)P;+1  +/P;-1 
'  2/  +  1 

equation  (58)  becomes 


(58) 


(59) 


(60) 


Since  this  is  an  infinite  series,  we  are  free  to  re-index  it  to  a  new  /,  and  then  take  two 
terms  in  the  series. 


-w^x  F-{ih+LLint} 

=  \v\Vr 


( Mi  ■  fP  ,  ILL A 

3  /o  1  5 


=  HVr(/0  COS  [0]+  /l/3) 

=  |v |  V  r  cos  [e  ]/0  +  |v|  V  r  •  7,  /3 
V  -V  rf  =  v-  Vr/0+vVr-7,/3 


(61) 


Where  we  made  use  of  the  definition  of  P, ,  and  the  term  f2  was  dropped  since  we  are 
considering  only  a  two-term  expansion. 

The  energy  flux  term  in  equation  (56)  is  tackled  in  a  similar  fashion,  although  we 
first  need  to  rewrite  the  expression  Vv .  Assume  we  align  a  hypothetical  z  axis  along  the 


same  direction  as  a ,  then  only  the  z  component  of  Vv  /  will  contribute  to  a  ■  Vv/  .  For 
that  case  we  write. 
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V  =^~ 
v  dv, 

vz  =vcos[0] 

y  _  dv  d  +  8cos[0]  d 


and  use  further  relations  to  recast  the  term 


dvz  dv  dvz  8cos[0] 
dv 


dv,  ’ 


(62) 


2  2  2  2 

VZ  =v  X+V  y+Vz 

|^-=^=oos[e] 
dv,  V 


We  also  recast 


8cos[0j 

dv , 


dv,  dv,  V  V  V  V  V 


then  combine  both  results  with  (62)  to  give 


a  mi  d  sin2[0]  9 

V=cos[0]— +- 


3v  v  8cos[0] 


(63) 


(64) 


(65) 


Continuing  the  assumption  of  alignment  along  a ,  and  using  the  results  above,  we  are 
able  to  write  the  energy  flux  term  in  (56)  as 


a-V  J  =d^—  =  a\ 

dv. 


cos[0]^  + 


df  ,  sin2[0]  df 


dv  v  8cos[0] 
Now  make  use  of  fundamental  identities  for  Legendre  polynomials 

cos[0]/; 

2/  +  1 

sin  + 

8cos[0]  21  +  1 


(66) 


(67) 

(68) 


Expanding  (66)  in  spherical  harmonics  and  using  the  identities  above  we  obtain, 
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a-Vv/  =a 


( a 


sin2[0]  3 


-cost^MlH  ^  3cos[ffJ1 


X/,f,[cos[S]] 


ravl 


f(l+DPM+lpM) 

|y/,py  +  l)W-l-^i))] 

21  +  1  J 

1  v{  21  +  1  [ 

(69) 


Since  the  summations  are  infinite,  we  are  again  free  to  re-index  them,  and  solve  for  the 
zeroth  and  first  order  terms. 


a-Vvf  =  a 


^ia/,-1  l  ,  a/j+1  /  +  !  ,  fM  (/  +  !)(/ +  2)  /w  Z(Z-1) 


=  « 


dv  2/-1  3v  21  +  3  v  21  +  3 


'JA+lA+*Amsm 

3  av  3  v  av 

'  i  a(vV,)+^.cos[ei 


v  21-1 


3v2  5v  dv 


Recalling  that  /,  is  a  vector,  the  final  result  for  the  energy  flux  is, 


(70) 


n  V  f-  -  d(v^~ 7i).^-va/o 
v/"3v2  av  v  av 


(71) 


The  final  term  in  (56)  is  the  collision  integral.  As  first  shown  by  Holstein  [10,  pp 
368-372],  this  has  two  components:  the  contribution  from  elastic  collisions  and  those 
from  inelastic  collisions.  Following  his  derivation  we  arrive  at  an  equation  for  each 
component. 


%~\  =  ~Nv  cos[0]  /,  Qe,as  [v] + ^-^-^-  (v  4/0  Qe,as  [v]) 

at  plastic  V  M  dv 

' collisions 


$■) 


.  =  2>,v 


fv’  \ 


/oK]1  Yi 


VV  ) 


G,^[v;]-/o[v]Qrw[v]-/1[v]cos[0]fir“[v] 


>  inelas 


inelas  r 


V nelasth 
'collisio, 

where  Q  refers  to  the  cross  sections  for  either  elastic  or  inelastic  collisions,  and  l  to  a 
specific  inelastic  process.  TV  is  the  colliding  partner  number  density,  v  the  velocity  for 


(72) 


(73) 
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which  the  distribution  function  is  being  solved,  i.e.  f0[r,v] ,  and  v't  a  different  velocity  at 
which  a  specific  inelastic  collision  /  is  taking  place. 

All  of  the  terms  are  now  combined  by  substituting  (61),  (71),  (72),  and  (73)  into 
(56).  After  moving  to  cylindrical  coordinates,  and  expressing  the  dot  product  using 
cos[6],  the  Boltzmann  equation  in  steady  state  becomes 

veos[fl]^  +  ^-j-(rfr)+_12  -d  (v2a  /i)  +  flCOs[e]|^-=  (74) 
dr  3  r  dr  3v  dv  dv 

-Nvcosmi  eeto[v]+4~T-(v4/0  Qe,asM) 

vMdv 

(  fv'  V  ^ 

+5>|V  /«[<]  -  2rte[v;]-/0[v]Grto[v]- u^cosmQ^iv] 

‘  [  \v ) 


We  can  split  the  above  equation  into  three  parts  by  first  integrating  over  6  to 
remove  the  cos[6\  terms;  this  gives  (75).  Next  we  multiply  (74)  by  cos[6]  and  again 
integrate  over  6;  this  removes  those  terms  without  a  cos[B\.  The  resulting  equation  is 
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The  next  step  is  to  change  variables  from  velocity,  v,  to  kinetic  energy,  u.  At  this 


€  -* 

time  we  also  include  the  form  of  the  acceleration  given  by  a  =  — -E .  With  these 

m 


changes,  the  three  equations  (75)-(77)  become, 

^{rfr)-e^[^Ur]fr  +  KfJ-^-[G[r,u]U]+uH[r,u]fQ=S0[r,u,f0}  (78) 

3r  or  3  du  du 

~ e0Er[r]^+K[r,u\fr  =0  (79) 

dr  du 

-eflEz^+K[r,u]fz=  0  (80) 

du 

Where  the  new  terms  below  follow  the  formalism  of  Uhrlandt  and  Winkler  [7],  and  are 
fully  described  in  section  3.1. 


G[r,M]  =  Y2^-M2A,ef[M]  (81) 

k  Mk 

H[r,u]^J^Nk  &7M  +  X2X  Qti  QtM  (82) 

k  l  k  l  k 

K[r,u]^Nk  Ql'[u]  +  H[r,u ]  (83) 

k 

5o^«./o]sZX «;  A7  /0[^«;i  (84) 

*  i 


In  writing  (81)-(84)  we  have  anticipated  the  need  to  consider  more  than  one 
neutral  species,  which  leads  to  the  summation  over  k.  The  S0  term  shown  is  simplified 

from  the  one  in  equation  (17).  Where  that  equation  is  specific  about  the  types  of  inelastic 
processes  considered;  equation  (84)  is  more  general. 

The  next  step  is  another  change  of  variables,  this  time  from  kinetic  energy,  u  to 
the  total  energy  e[r].  The  total  energy  is  a  function  of  r  because  of  the  radial  electric 
field.  The  relationship  is 
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£[r]  =  u-e0V[r] 

(85) 

r 

V[r]  =  -\Er[r]dr 

(86) 

0 


With  this  transformation,  the  distributions  are  functions  of  e,  and  we  adopt  the  tilde  to 
differentiate  the  two  forms. 


fx[r,u[r,£]]  =  fx[r,£] 

—f  =— f 
duJx  de  x 


(87) 

(88) 
(89) 


where  x  can  refer  to  0,  r,  or  z. 

With  the  transformation  to  total  energy,  equations  (79)  and  (80)  are  solved  for 


fr  and  fz  respectively  to  give 


K[r,e]  de 

We  then  use  (90)  and  (91)  to  eliminate  fr  and  fz  from  (78).  The  result  is  the  elliptic 
partial  differential  equation, 


(90) 

(91) 


ru 


-/o)+T— ( 


a  MeaEz )2  a 


r  dr  3AT[r,fr]  3r  de  3 K[r,e]  de  de 


fo)  +  ^(Glr,e]fo)-uH[r,e]fo+So[r,eJ0]  =  O  (92) 


This  is  the  equation  used  in  the  strict  solution. 
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Appendix  C  Derivation  of  Upper  Boundary  Condition 


This  appendix  derives  an  approximate  EEDF  based  on  a  limiting  form  of  equation 
(10),  valid  at  high  energies.  The  resulting  solution  is  used  as  an  upper  boundary 

condition  for  /0  !/,£„,] ;  this  boundary  condition  is  discussed  further  in  section  3.1.3  page 


26. 


The  derivation  begins  with  the  strict  kinetic  CBE, 

(93) 

At  the  highest  energies  the  system  is  approximated  by  a  balance  between  the  field  driven 
flux  and  the  losses  from  inelastic  collisions.  The  equation  takes  the  form 

(94) 

de  3 K[r,e]de  °  0 

We  assume  an  EEDF  of  the  form  f0=oc  exp[-/3  u] ,  substitute  this  into  (94),  and  expand 
the  result. 

ZS  ..  O  ,  .TTV  „ 

(95) 


d  (  u  R  ~  3uH[r,e]  ~ 
\  _-A  P)J o'~  ,  x2  Jo 


/o 


^  p2u 


de  K[r,e] 

a f 


0 e.Ezy 


-P- 


isT[r,e]  de 


W 


K[r,e] 


To  continue,  we  must  assume  a  functional  form  for  the  inelastic  cross  section.  In  the  case 
of  neon,  assume  a  weak  dependence,  Q[u]  =  Q0u'n .  Then 


a 

f  u  ] 

_  a 

r  U  ^ 

_  a 

{  1 

2  „  2  (96) 

de 

[K[r,e]  J 

”  de 

[a0[u]  J 

de 

3umNQ0  3  K[r,e] 

After  substituting  back  into  (95),  solve  for  (3. 
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f^u  2/3 
K[r,e]  3  K[r,e] 


3 uH[r,e]  ~ 

17aF/o 


M 

K[r,e] 


[fiu-2/3]= 


7o  fu  _ 

*[r,£] 


^2=^4rr^r’e^r’£] 

(«A) 


(97) 


where  the  term  2/3  is  dropped  by  assuming  it  is  small  in  comparison  to  j3u.  A  quick 

check  of  this  assumption  is  appropriate.  For  typical  values  in  a  Neon  discharge, 

Neutral  density  ~  2  •  1016  cm"3  Qelastic  ~  3  •  10"16  cm2 

Ez  ~  2V/cm  Qinelastic~1.4-10-17cm2 

“max  ~  30eV 

/3  ~  0.4 ,  and  thus  j3  u  ~  12.0 ,  which  justifies  the  earlier  assumption. 

Based  on  the  earlier  assumed  form  for  the  distribution,  f0=oc  exp [~/3  u] ,  we  can 
relate  f0[r,jmax]  to  f0[r,jmx-l],  where y  signifies  an  energy  bin  counter.  To  simplify  the 
expression  assume  (3[u]  does  not  vary  between  jmax-l  and  jmax,  then 

7o  [r,  jmm  ]  =  A  7,  [r,  ymax  - 1]  (98) 

A  =  exp[-/j  (nymax  -Mymax_,)] 

VMJ 
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Appendix  D  Boltzmann  Equation  in  Finite  Differenced  Form 


Finite  differencing  begins  with  the  analytic  CBE  derived  in  Appendix  B, 


For  the  computer  implementation,  we  change  notation  slightly.  The  two  functions  K[r,e ] 
and  H[r,£]  are  given  descriptive  names,  mtr[r,e]  and  inelas[r,e]  respectively  (see 
discussion  in  section  3.1.1).  Furthermore,  let  /„—»/,  to  ease  the  notation.  With  these 
changes,  rewrite  (99)  as, 

-^-(rH[r,£]^-f)+^-(D[r,e]~f)  +  ^-(G[r,£]f)-uinelaS[r,£]f  +  S0[r,e,f]  =  0  (10 

r  dr  dr  d£  d£  d£  0) 

where 


A  11 

H[r,£\~ 

3  mtr[r,£] 

(101) 

u(e  E  )2 

D[r,g]^/  0  *} 

3  mtr[r,£] 

(102) 

Now  tackle  the  terms  in  (100)  individually.  First  the  radial  flux  term, 

1  ^  0  ^ 

(rtf[r,£]— /0) .  Using  the  differencing  equations  (36)-(40),  we  write  this  as 
r  dr  or 

~(r»[r,£)|:/)srw/WJ-(rw  +  S1J)/1J+Su/HJ  (103) 

T,  J  =  (<*,  (2 + a,  k,  Hmj  +  (l -af  )rfi, ,,  +  ) 

S,,js— L  +  (l  -  A2  )rfiLJ  +  A  (2  +  A  k,  H,_, j ) 

'i 

Next  the  energy  terms,  using  modified  forms  of  (36)-(40): 
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(104) 


and 


^(Dtr,£]^/)SWu/w  -<TW  +Y,jf„H 


W,J  -  c)  (fj  (2 + Yl  )Dut,  +  (1  -  y)  )d,j  +  d,m  ) 
Kj*c;(DIJtl +(1  -«;K, 


— (G[r,£]/0)  3  K,Jut,+(Nu -L^lj 


k,., 

'  i' <v 


AT.  =C  8  G  . 

«.j  7  7  1.7 


(105) 


As  the  last  two  terms  of  (100)  don’t  involve  any  derivatives,  they  are  evaluated 
simply  at  each  i,j.  The  equations  given  above  form  the  complete  set  for  all  interior 
points  on  the  grid.  At  the  boundaries  we  must  use  forward  or  backward  differencing,  and 
consider  the  boundary  conditions. 

Start  first  with  the  r  =  0  boundary.  Note  that  the  only  term  that  needs 
modification  is  the  radial  flux  term,  since  as  (103)  is  written,  it  tries  to  access  a  point 
ft_X  j  outside  of  the  grid.  To  fix  this  we  first  take  the  limit  of  the  radial  flux  as  r  — >  0, 


1  3  1  d  j, . 

~—(rH[r,e]  —  f) 
r  dr  dr 


=  2  H[r,e] 


a2/ 


r  — >  0 


dr 1 


r  =  0 


(106) 


Since  we  are  on  the  edge  of  the  grid,  invent  a  hypothetical  grid  point  at  i  =  -1  and 


assume  uniform  grid  spacing,  then  d2/  becomes 

dr 2 


a2/  ~  fi-u  -2fi,j+fM,j  (107) 

dr 2  A  r 

Use  the  r  =  0  boundary  condition 
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dr 


r= 0 


=  0  ^IldZlihL 

2rx 


2H[r,e] 


ki  =  f-hj 

a2/ 1  _4/f0J(4,-/0,;.) 


3r2 


(108) 


r=0 


At  the  next  boundary,  u  =  0,  directly  enforce  the  boundary  condition  given  in 
(19).  The  condition  is 


8<p[r]  3/ 


«= o 


3r  3r 


=  0 


(109) 


u=0 


where  <p[r]  =  -V[r].  Since  both  derivatives  involve  the  boundary,  neither  can  be  treated 
with  central  differencing.  Instead,  employ  backward  differencing  for  the  radial 
derivative,  and  forward  for  the  energy  derivative.  As  discussed  in  section  3.2.3,  we  have 
made  a  decision  to  limit  our  radial  increments  to  +/-  1  around  these  solution  points,  in 
order  to  speed  the  solution.  This  forces  us  to  use  a  two-point  derivative  in  the  radial 
direction,  while  for  the  energy  derivative  we  are  free  to  use  equation  (37).  The  resulting 
differenced  equations  are 

o  =  -xi,  kj  +  R2,  A,.+1  -  R3,  fiJ+2  +  RA.t  /M,,  (110) 

( d(p 


Rl  = 


dr 


4-^+ e;cm(sM+2) 
An-n-il 


R2t=ElCM(y„+SM+2) 

R3i  —  Ci+iYi+i 

1 


,•  (n-n-i) 
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There  is  a  special  point  at  the  origin  where  both  the  r  =  0  and  u  =  0  boundary 


Z\f 

conditions  hold.  At  this  point  —  =  0  and  (110)  reduces  to 

dr 


0  =  -0O  +  Q20  /0>1  -  Q\  /0;2  (111) 

do ft +2) 

Q20  +2) 

Q30  =  £z  Q/j 

At  the  wall  boundary,  similar  to  the  r  =  0  boundary,  only  the  radial  equation 
needs  alteration.  The  technique  is  to  substitute  for  fi+l  j  in  equation  (103)  with  a  discrete 
version  of  whatever  wall  loss  condition  is  used.  In  the  case  of  the  wall  loss  function,  this 
leads  to 


~(rH[r.«]|:/)sr,JFu -<JtJ  +  Su)/„  +(T(J  +SU)/HJ 
s-2(r,-r,_l)  mtr[ri,£i j\A  exp [aufj] 

where  A  and  a  are  input  parameters  as  shown  in  equation  (21). 

For  the  wall  loss  cone,  the  radial  flux  becomes 

i  ±(rHlr,e]^-f)  =  -(,Tl.lX,J+TliJ  +S„)/U  +(T„  +S, 


xij-2(ri  ~ri- i) 


3(1-*) 

2(1+*) 


(112) 


When  j  =  NJ- 1,  we  substitute  into  (104)  and  (105)  for/-J+;.  Using  equation  (20), 
the  energy  flux  terms  become 


* 
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(113) 


Be 


(G[r,e]f0)  =  (Kt  jAA+N;j  —LUJ)fitj  +Mijfij_l 


where 


AA  =  exp  [-6  (u  —  u  ,Y| 

wvrL  h'  K^j  max  y  max-1 


/3  = 


//ifr  [r,  £max  ]  me/<w[r,£maJ 


This  completes  the  set  of  differenced  equations. 


(114) 


(115) 
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